AD-A192  212 


AFGL-TR-87-0250 


«TH»  FflT  COPv 


Forward  and  Inverse  Modeling  of  Near-Field 
Seismic  Waveforms  From  Underground  Nuclear 
Explosions  for  Effective  Source  Functions 
and  Structure  Parameters 


L .  J .  Bu  rd  i  ck 
J.S.  Barker 


Woodward-Clyde  Consultants 
566  El  Dorado  Street 
Suite  100 

Pasadena,  CA  91101 


5  April  1987 


DTIC 

ELECTE 
FEB  0  3  1988 


Final  Report 

5  September  1986-5  March  1987 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


AIR  FORCE  GEOPHYSICS  LABORATORY 
AIR  FORCE  SYSTEMS  COMMAND 
UNITED  STATES  AIR  FORCE 

HANSCOM  AIR  FORCE  BASE,  MASSACHUSETTS  01731-5000 


88  1 


'i  6  U50 


This  technical  report  has  been  reviewed  and  is  approved  for  publication." 


HENRY  A.  06SING 


4. 


Chief,  Solid  Earth  Geophysics  Branch 


FOR  THE  COMMANDER 


JONALD  HV  ECKHARDT 
Director 

Earth  Sciences  Division 


This  report  has  been  reviewed  by  the  ESD  Pub? Ic  Affairs  Office  (PA)  and  is 
releasable  to  the  National  Technical  Information  Service  (NTIS). 


Qualified  requestors  may  obtain  additional  copies  from  the  Defense  Technical 
Information  Center.  All  others  should  apply  to  the  National  Technical 
Information  Service. 


If  your  address  has  changed,  or  if  you  wish  to  be  removed  from  the  mailing 
List,  or  if  the  addressee  is  no  longer  employed  by  your  organization,  please 
notify  AFGL/DAA,  Hanscon  AFB,  MA  01731.  This  will  assist  us  in  maintaining 
a  current  mailing  list. 


Do  not  return  copies  of  this  report  unless  contractual  obligations  or  notices 
on  a  specific  document  requires  that  it  be  returned. 


UNCLASSIFIED 


S£CU«'TV  Classification  Of  This  PAGE 


REPORT  DOCUMENTATION  PAGE 

U  REPORT  se  CuRiTV  CLASSIFICATION 

Unclassified 

lb  RESTRICTIVE  MARKINGS 

2*  SECURITY  CLASSIFICATION  AUTHORITY 

3  DISTRIBUTION.  AVAILABILITY  OF  REPORT 

Approved  for  public  release;  distribution 

2b  DEC  LASS  F  .CATION 'DOWNGRADING  SChEOULE 

uni imited. 

4  PERFORMING  ORGANIZATION  REPORT  NUM8ERIS) 

.  WCCP-R-87-91 

5.  MONITORING  ORGANIZATION  REPORT  NUM8ERIS) 

AFGL-TR-87-0250 

6«  NAME  of  PERFORMING  ORGANIZATION 

Woodward-Clyde  Consultants 

6b.  OFF  ICE  SYMBOL 
(If  applicable  > 

7*.  NAME  OF  MONITORING  ORGANIZATION 

Air  Force  Geophysics  Laboratpry 

6c  AOORESS  . City  State  and  ZIP  Code  l 

7b  ADDRESS  (City.  State  and  ZIP  Code t 

566  El  Dorado  Street,  Suite  I 
Pasadena,  CA  91101 

O 

O 

Hanscom  Air  Force  Base 

Massachusetts  01731 

Ba  NAME  OF  F  UN  O  t  NG /SPONSORING 
ORGANIZATION 

DARPA 

0b  OFF  ICE  SYMBOL 
<lf  applicable ) 

9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

F19628-85-C-0036 

8c  ADDRESS  City  Stale  and  ZIP  Code) 

1400  Wilson  Boulevard 
Ar 1 i ngton ,  VA  22209 


11  T  TWE  -Inciude  Security  Classification) 

See  block  16;  unclassified 


12  PERSONAL  AjTnORlSi 

L.  J.  Burdick  and  J.  S.  Barker 


10  SOURCE  OF  FUNDING  NOS 


1  3a  T  v  p £  OF  REPORT 

Final  Report 


16  SuPPLEMENTAP  Y  NOTATION 


13b  TIME  COVERED 


14  OATE  OF  REPORT  i  Yr  .  Mo  .  Day  I  I  15  PAGE  COUNT 


87/4/5 


Forward  and  Inverse  Modeling  of  Near-Field  Seismic  Waveforms  from  Underground  Nuclear 
Explosions  for  Effective  Source  Functions  and  Structure  Parameters _ _ 


COS  AT  i  CODES  [  18  SUBJECT  TERMS  (Continue  on  nri«rt«  if  neceuary  and  identify  by  block  number) 

Inversion  of  explosion  seismic  waveforms,  near-field  obser¬ 
vations,  effective  source  functions,  strain,  non-linear 
material  response  _ _ _ 


A  AS  T  PACT  'Continue  on  neue^te  if  neceteary  and  identify  by  bloc*  number/  ]  ACC09SiOn  FOT 


19  ABSTRACT  'Continue  on  reverie  if  neceuary  and  identify  by  block  numberi 

See  reverse  side 


nt;s  gfasj. 

DTI  ?  T A >< 
lUittrrjtftau'Cd 
J ; .  -.i  l  ;  f  1  n :  '  1 


■  List  ci but i on/ 

AvaiU'.bl Jlty  Cedes 
Avail  aiid/or 


V  * 


Da  st.  ( 

|  Cpec: 

!  1 
! _ L 

20  OiS'P  fluT'ON  A.AaABiL  Tv  OF  ABSTRACT 

-M- 

jnCw ass  (  e  ::  us.  w  Tf  :  same  as  «pt  S:  otic  users 


?2»  SAME  SF  PESPGNS  b^E  iNCviDUAt 

Jar es  Lewkowicz 


22b  telephone  number 

i Include  .4  < to  Code 

(617)  377-3028 


22c  OFF  'CE  S  vmBOl 

AFGL/LWH 


UNCLASSIFIED _ 

.  SECURITY  CLASSIFICATION  Of  THIS  PAGE _ _ 

\ 

It  is  well  established  that  near  field  records  of  nuclear  explosions  can  be  analyzed  to 
obtain  detailed  information  about  the  seismic  source  function  and  its  dependence  on  yield. 
This  information  is  generally  formulated  in  terms  of  parameterized  models  for  the  RDP  and 
for  the  test  site  crustal  structure.  In  this  study , .review  the  results  of  forward  model¬ 
ing  studies  to  obtain  the  source  and  structure  parameters  for  Pahute  Mesa.  These  models 
fit  the  observed  near  field  records  well,  but  there  is  some  question  as  to  how  errors  in 
crustal  structure  might  affect  seismic  source  parameters.  Furthermore,  there  is  a  general 
need  to  be  able  to  develop  source-structure  models  in  a  consistent,  unbiased  fashion.  To 
address  these  issues,  we  have  developed  a  simultaneous  inversion  for  source  and  structure 
parameters.  In  previous  reports,  we  discussed  the  development  of  the  method  and  applied 
it  to  Pahute  Mesa  data.  In  this  report,  we  present  an  application  of  an  improved 
technique  for  inverting  for  parameters  in  these  types  of  problems  known  as  the  jumping 
method.-  The  primary  advantage  of  this  method  is  that  parameter  constraints  and  error 
analyses  are  made  on  the  parameters  of  the  models  rather  than  on  the  changes  in  the 
parameters.  We  show  results  of  application  of  the  method  to  artificial  and  real  near¬ 
field  data  sets.  Both  tne  utility  and  limitations  of  the  method  have  been  eluciadated. 
Tradeoffs  between  parameters  are  quantified  by  the  error  analysis  inherent  in  any 
inversion  method. >  A  second  problem  with  analysis  of  near-field  records  from  explosions 
is  that  there  is  some  question  as  to  whether  crustal  materials  respond  in  a  linear 
anelastic  or  anelastic  fashion  or  whether  they  have  significant  nonlinear  response 
because  of  the  high  shear  strain  levels.  Near-field  data  sets  from  several  nuclear 
explosions  and  an  earthquake  are  examined  to  address  this  question.  It  is  shown  that  if 
nonlinear  naterial  response  occurs  it  does  not  have  a  large  enough  effect  to  have  a 
significant  effect  on  data  interpretation.  . 
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INTRODUCTION 


Perhaps  the  most  crucial  step  in  the  determination  of  effective  source 
functions  and  the  estimation  of  yields  from  near-field  recordings  of 
underground  nuclear  explosions  is  the  development  of  an  accurate  velocity 
structure  model.  For  ranges  outside  of  the  spall  zone,  the  first  peak  in 
observed  vertical  velocity  records  is  often  the  result  of  the  interference 
between  the  direct,  upgoing  wave  and  energy  that  departs  downward  from  the 
source  but  is  turned  by  the  velocity  gradient  in  the  upper  crust.  Thus  the 
amplitude  and  width  of  the  source  pulse  inferred  from  the  observed  waveform  is 
quite  sensitive  to  the  velocity  gradient,  particularly  near  the  source  depth. 
Previous  studies  have  obtained  detailed  velocity  structure  models  through 
trial-and-error  waveform  modeling  (e.g.,  Helmberger  and  Hadley,  1981; 
Hartzell,  et  al.,  1983;  Burdick,  et  al.,  1984).  We  have  used  the  model  of 
Hartzell,  et  al  (1983)  in  a  forward  modeling  study  for  the  source  parameters 
of  four  Pahute  Mesa  explosions.  The  resulting  synthetic  seismograms  provide 
an  admirable  fit  to  the  observed  waveforms.  There  is  some  concern,  however, 
regarding  the  adequacy  of  the  structure  model  for  a  variety  of  locations 
within  Pahute  Mesa.  In  particular,  we  need  to  consider  how  potential  errors 
in  the  structure  model  effect  the  values  determined  for  the  source  parameters. 
An  additional  problem  is  that  in  determining  source  parameters  from  explosions 
in  a  completely  different  structure  (e.g.,  Yucca  Flats,  Amchitka,  or  one  of 
the  Soviet  test  sites),  a  new  velocity  structure  model  must  first  be 
determined . 


In  order  to  address  these  issues,  we  have  developed  a  simultaneous 
inversion  method  which  solves  for  the  parameters  of  both  the  source  and 
velocity  structure  models.  A  report  on  the  initial  development  of  the  method 
was  given  in  Barker,  et  al .  (1985b)  and  the  application  of  this  method  to 
Pahute  Mesa  data  was  presented  in  Barker,  et  al .  (1986).  We  have  improved  the 
inversion  technique  by  incorporating  the  jumping  method  of  Shaw  and  Orcutt 
(1985).  The  primary  improvement  in  this  technique  is  that  parameter 
constraints  and  error  analyses  are  made  on  the  parameters  of  the  model,  rather 
than  on  the  changes  in  the  parameters.  We  present  a  test  inversion  of  a 
synthetic  data  set.  The  results  are  quite  promising,  but  also  illustrate  the 
limitations  in  the  resolution  of  deeper  structures  when  limited  ranges  of 
observations  are  considered. 

The  inversions  of  waveforms  for  same  four  Pahute  Mesa  explosions  are 
presented  and  compared  with  the  forward  modeling  results.  In  general,  the 
synthetic  waveforms  from  the  inversions  provide  improvements  over  the  forward 
modeling  in  amplitude,  arrival  time,  and  in  some  cases,  in  the  detailed  shapes 
of  the  waveforms.  For  BOXCAR,  the  structure  model  obtained  is  not 
significantly  different  than  that  of  Hartzell,  et  al .  (1983).  This  is  not 
particularly  surprising,  since  these  data  were  modeled  in  the  development  of 
the  Hartzell,  et  al .  (1983)  structure.  For  the  other  events  modeled,  there  is 
very  little  resolution  of  structure  parameters  deeper  than  about  3  km.  Above 
this,  however,  we  will  show  that  the  velocity  gradient  near  INLET  is  verv 
similar  to  that  near  BOXCAR,  but  the  gradient  is  much  greater  near  MAST.  The 
results  from  SCOTCH  suggest  that  the  water  table  is  somewhat  deeper  than  the 


constrained  value  'used,  but  that  the  gradient  below  this  is  similar  to  that  of 
BOXCAR.  The  variation  in  velocity  gradient  for  MAST  may  be  related  to  its 
location  on  the  northeast  edge  of  Pahute  Mesa. 

From  the  parameters  of  the  effective  source  functions  we  may  develop 
scaling  relations  and  estimate  yields  for  the  events  studied.  Barker,  et  al . 
(1985a)  estimated  yields  from  the  forward  modeling  results.  In  this  report, 
we  use  two  new  relations  which  are  constrained  to  the  announced  yields  of 
BOXCAR  and  SCOTCH.  Yields  are  estimated  for  all  four  of  the  events  for  which 
the  simultaneous  inversion  has  been  performed. 

While  the  resolution  of  the  inversion  results  depends  on  the  range 
distribution  of  the  observations,  and  the  results  are  limited  by  noise  or  the 
effects  of  lateral  structural  variations  in  the  observed  data,  the  application 
to  Pahute  Mesa  waveforms  has  demonstrated  the  usefulness  of  the  method. 
Trade-offs  between  parameters,  particularly  between  structure  and  source 
parameters,  are  quantified  by  the  error  analysis  inherent  in  the  inversion 
method.  Finally,  the  inversion  may  be  applied  to  near-field  waveforms  from 
any  site,  without  any  a  priori  information  on  the  crustal  velocity  structure. 
This  makes  such  an  inversion  valuable  as  one  of  many  on-site  verification 
techniques,  enabling  the  calibration  of  test  sites  for  treaty  monitoring 
purposes . 

NEAR-FIELD  DATA  SET 

The  near-field  data  discussed  in  this  study  consist  of  three -component 
surface  recordings  of  acceleration  and  velocity  digitized  by  Sandia  National 


Laboratories  and  reported  by  Perret  (1976).  The  events  modeled  are  listed  in 
Table  1  and  include  BOXCAR,  MAST.  INLET  and  SCOTCH.  The  distribution  of 
stations  for  each  event  is  given  in  Table  2  and  displayed  in  Figure  1. 


Table 

1  -  Pahute 

Mesa  Events 

Modeled 

Event 

Date 

Depth 

Yield 

(m) 

(kt) 

BOXCAR 

4/26/68 

1160 

13001 

MAST 

6/19/75 

912 

5202 

INLET 

11/20/75 

817 

5002 

SCOTCH 

5/23/67 

970 

1551 

1  announced  yield  (Marshall,  et  al . ,  1979) 

2  estimated  yield  (Dahlman  and  Israelson,  1977) 
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Table  2  -  Station  Distribution  and  Arrival  Times 


BOXCAR 

Range 

Arrival  time1 

Residual3 

Window  start4 

Window  end4 

(km) 

(sec) 

(sec) 

(sec) 

(sec) 

S  - 12 

3.81 

1.26 

0.09 

1.25 

2.05 

S-16* 

4.87 

1.51 

0.16 

1.55 

2.35 

S-24* 

7.27 

2 . 10 

0.16 

1.85 

3.20 

S  -  34* 

10.37 

2.92 

-0.03 

2  .  70 

3.95 

S  -  74* 

22.47 

4.92 

0.32 

4.50 

6.00 

MAST 

Range 

Arrival  time2 

Residual3 

Window  start4 

Window  end4 

(km) 

(sec) 

(  sec ) 

(sec ) 

( sec ) 

S  -  5* 

3.65 

1.15 

0.14 

1.00 

1.90 

S  -  6 

5.47 

1.56 

0.26 

1.40 

2.35 

S  -  7 

7.30 

2.03 

0.29 

1.85 

2.80 

INLET 

Range 

Arrival  time2 

Residual3 

Window  start4 

Window  end4 

(km) 

(sec) 

(sec) 

(sec) 

( sec ) 

S  -  5 

1.63 

0.59 

0.09 

0.40 

1.25 

S-6* 

3.27 

1.03 

0.15 

0.90 

1  .  75 

S  -  7* 

6.53 

1.99 

0.16 

1.85 

2  .  70 

SCOTCH 

Range 

Arrival  time1 

Residual3 

Window  start4 

Window  end4 

(km) 

(sec) 

(sec) 

(sec) 

( sec ) 

S-3A* 

4.13 

1.38 

0.06 

1  .  20 

2 . 00 

S  -  4* 

6.06 

1.89 

0.12 

1  .  75 

2  .  70 

*  Station  used  by  Hartzell,  et  al.  (1983)  to  determine  the  structure  mode 
Station  located  outside  Silent  Canyon  caldera 

*  Arrival  times  from  Perret  (1976) 

2  Measured  arrival  times 

3  Time  predicted  by  Hartzell,  et  al .  (1983)  model  -  observed  time 

*  Time  windows  used  in  the  waveform  inversion 


PAHUTE  MESA  VELOCITY  STRUCTURE  MODEL 


The  velocity  structure  model  for  Pahute  Mesa  was  obtained  by  modeling  the 
arrival  times  and  waveforms  of  the  surface  velocity  records  indicated  in  Table 
2.  The  model  was  previously  presented  by  Hartzell,  et  al .  (1983)  and  is 
listed  in  Table  3.  A  plot  of  the  model  as  a  function  of  depth  is  shown  in 
Figure  2,  along  with  a  comparison  with  other  models  for  Pahute  Mesa.  The  P- 
and  S-wave  models  of  Helmberger  and  Hadley  (1981)  are  shown  as  dotted  lines 
and  the  P-wave  model  of  Hamilton  and  Healy  (1969)  and  the  S-wave  model  of 
Carroll  (1966)  are  shown  as  dashed  lines. 


The  validity  of  the  model  may  be  evaluated  based  on  P-wave  travel  times 
and  the  fit  of  synthetic  waveforms  to  the  data.  The  S-wave  velocity  structure 
is  based  primarily  on  the  amplitude  and  arrival  time  of  the  long-period 
Rayleigh  wave.  The  ratio  of  vertical  to  radial  amplitudes  may  also  be  used  as 
a  criterion  since  it  is  particularly  sensitive  to  the  velocity  gradient  near 
the  surface.  However,  the  velocities  of  the  uppermost  layers  may  be  expected 
to  vary  laterally  across  Pahute  Mesa,  so  the  velocity  structure  model 
determined  must  be  considered  an  average,  not  intended  to  fit  variations  in 
the  amplitude  ratios. 


The  observed  P-wave  arrival  times  and  the  time  residuals  (predicted  - 
observed)  due  to  the  velocity  structure  model  are  listed  in  Table  2.  For  the 
nearer  receivers  the  residuals  are  all  less  than  0.1  s.  The  variation  in 
these  residuals  reflects  the  lateral  velocity  variations  in  the  uppermost 
layers.  Orly  the  receivers  at  ranges  greater  than  15  km  have  residuals  greater 
than  03  s  indicating  that  the  lower  layers  may  be  slightly  too  fast.  These 


Table  3  -  Pahute  Mesa  Velocity  Structure  Model* 


Layer 

P  Velocity 
( km/s ) 

S  Velocity 
(km/s ) 

Dens i ty 
(g/cm3) 

Thickness 

(km) 

1 

2  .  30 

1.35 

1.90 

0.36 

2 

2.80 

1.50 

2.00 

0.40 

3 

3.30 

1.50 

2.25 

0.70 

4 

4.00 

1.90 

2.30 

0.70 

5 

4.60 

2.00 

2.40 

0.75 

6 

5.30 

2.50 

2.50 

0.80 

7 

5  .  50 

2.95 

2.70 

2.25 

8 

6.10 

3.50 

3.00 

10.00 

9 

7.00 

4.00 

3.00 

10.00 

*  from  Hartzell,  et  al .  (1983) 
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receivers  are  also  located  outside  of  the  Silent  Canyon  Caldera.  Other  than 
two  rather  noisv  receivers  for  MAST,  all  receivers  within  10  km  have  residuals 
less  than  0.2  s,  indicating  that  the  model  fits  arrival  times  exceptionally 
we  1 1  . 


The  fits  of  the  synthetic  seismograms  to  the  observed  data  are  the 
subject  of  the  following  sections,  but  an  example  of  the  ability  of  the 
velocity  structure  model  to  account  for  extensive  portions  of  the  observed 
vertical  waveforms  is  shown  in  Figure  3.  Shown  are  the  vertical  velocity 
recordings  for  BOXCAR  at  ranges  of  7.26,  7.70  and  10.36  km,  along  with 

synthetic  seismograms  generated  using  the  wavenumber  integration  technique  of 
Yao  and  Harkrider  (1983).  Details  of  the  first  1. 5-2.0  s  after  the  P  arrival 
are  quite  well  modeled,  as  are  the  amplitude  and  arrival  time  of  the  Rayleigh 
wave  (indicated  by  an  arrow).  Later  P  arrivals  may  reflect  unmodeled 
phenomena  such  as  spall  or  tectonic  release.  The  observed  P  wave  at  10.36  km 
includes  frequencies  higher  than  those  computed  for  the  synthetics,  which 
accounts  for  the  amplitude  discrepancy  in  the  figure.  Figure  U  shows  a 
comparison  of  the  observed  waveforms  at  10.36  and  7.26  km  with  synthetic 
seismograms  computed  by  the  generalized  ray  method  of  Helmberger  and  Harkrider 
(1978).  While  these  synthetics  cannot  match  the  entire  waveform,  they  contain 
the  high  frequencies  required  to  model  the  first  several  arrivals.  These 
synthetics  will  be  used  in  subsequent  sections  for  determining  source 
parameters . 

Since  the  determination  of  source  parameters  will  depend  to  a  great 


extent  on  the  amplitude  and  width  of  the  first  pulse  of  the  observed  vertical 
velocity  waveform,  it  is  instructive  to  investigate  how  the  interaction  of 


Figure  3.  Comparison  of  the  vertical  velocity  records  at  distances  of  7.26, 
7.70,  and  10.36  km  for  BOXCAR  (top  traces)  with  synthetics  (bottom  traces)  for 
the  preferred  structure  in  Figure  2  using  the  full  wave  theory  code  KB. 
Amplitudes  in  cm/sec  are  peak  values  for  the  entire  record. 
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different  rays  through  the  structure  affects  the  waveshape.  Figure  5  shows 
the  decomposition  of  the  vertical  component  synthetics  into  various  ray  tvpes 
for  ranges  of  2-12  km  from  a  source  with  the  parameters  of  BOXCAR.  The  source 
is  located  at  a  depth  of  1.160  km.  For  each  range,  the  top  plot  consists  of 
the  upgoing  P  wave  only.  The  next  trace  is  the  sum  of  diving  P  waves  and 
consists  of  reflections  and  refractions  from  interfaces  below  the  source.  The 
third  trace  is  the  sum  of  pP  rays  that  reflect  from  the  free  surface  and  then 
reflect  or  refract  from  layer  interfaces.  Below  this  is  the  sum  of  P-to-S 
conversions  at  layer  interfaces  for  direct  rays,  or  at  reflection  points  for 
turning  rays.  Finally,  the  bottom  trace  for  each  range  consists  of  the  sum  of 
all  of  these  rays.  Each  trace  is  scaled  to  the  first  peak  amplitude  of  the 
sum  trace . 


At  2  km,  the  upgoing  P  wave  is  the  first  arrival  and  the  only  significant 
contribution  to  the  sum.  By  U  km,  the  diving  P  wave  is  beyond  the  critical 
angle  of  the  first  reflector  below  the  source,  its  amplitude  is  comparable  to 
that  of  the  direct  ray,  and  it  appears  impulsive  due  to  the  pos  t  -  c  r  i  t  ical 
reflection  phase  shift.  At  6  km,  the  diving  P  wave  is  the  first  arrival  and 
consists  of  two  post-critical  reflections,  direct  P  serves  only  to  broaden  the 
pulse,  and  pP  is  beginning  to  emerge.  At  8  km,  diving  P  is  the  earliest  and 
largest  phase,  direct  P  destructively  interferes  with  the  backswing  and  pP  is 
an  impulsive  late  arrival.  At  10  km  and  beyond,  individual  diving  rays  begin 
to  separate  in  time  and  interfere  with  pP,  which  now  consists  of  multiple 
pos t - cr i t ical  rays.  For  this  velocity  structure  model,  at  all  ranges,  P-to-S 
conversions  contribute  only  negligibly  to  the  vertical  synthetics. 
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To  summarize  the  effects  of  the  velocity  structure  model  on  the  synthetic 
waveforms,  the  upgoing  P  wave  alone  is  sufficient  only  for  ranges  less  than 
about  4  km.  In  practice,  most  vertical  waveforms  for  these  near  ranges  would 
be  truncated  by  the  onset  of  spall  immediately  after  the  first  peak.  For 
ranges  of  about  4-8  km,  the  interference  of  upgoing  and  diving  P  ravs  can 
greatly  influence  the  amplitude  and  pulse  width  of  the  first  peak.  Similarly 
for  ranges  greater  than  8  km,  the  first  peak  is  sensitive  to  the  interference 
of  various  post  -  critical  diving  rays,  and  the  scale  of  discretization  of  the 
velocity  gradient  becomes  important.  pP  contributes  complications  later  in 
the  P  waveform  for  ranges  beyond  about  8  km,  but  S-to-P  conversions  are  not 
significant  in  the  vertical  synthetics  at  any  range. 

An  implicit  assumption  in  this  modeling  approach  is  that  the  velocity 
structure  may  be  parameterized  as  a  stack  of  horizontal,  homogeneous  layers. 
The  effects  of  discretizing  the  vertical  velocity  gradient  have  been  discussed 
elsewhere  (e.g.  see  Burdick,  et  al.  1982),  but  in  general,  as  long  as  layer 
thicknesses  are  less  than  the  predominant  wavelength  of  the  observat ions ,  the 
discretization  will  not  adversely  affect  the  results.  While  lateral 
variations  in  the  velocity  of  the  layers  above  the  source  or  in  the  thickness 
of  the  top  layer  due  to  changes  in  elevation  may  be  noticeable  in  arrival 
times  and  in  P-to-S  conversions  on  the  radial  component,  these  variations  do 
not  appear  to  have  a  significant  effect  on  the  vertical  waveforms.  On  the 
other  hand,  the  plane-layer  assumption  tends  to  break  down  as  range  increases 
beyond  about  10  km,  due  to  the  sensitive  interaction  of  head  waves  and 
post-critical  reflections  from  deeper  layers.  Finally,  some  of  the  receivers 
modeled  in  this  study  lie  outside  the  boundary  of  Silent  Canyon  Caldera,  the 
geologic  feature  that  dominates  the  subsurface  structure  of  Pahute  Mesa.  For 
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WAVEFORM  MODELING  FOR  SOURCE  PARAMETERS 


The  vertical  particle  velocity  due  to  an  explosion  source  may  be  written 
using  first  order  asymptotic  generalized  ray  theory  as  (see  for  example, 
Helmberger  and  Harkrider,  1978) 


I  2d  1  r^s 

u  (R.t)=  W(ty  ,  --  T*  ^  R,(P) 
n  \  Rdt\,  t 


where  *  denotes  convolution,  ^ ( t )  is  the  second  time  derivative  of  the 
reduced  displacement  potential  (RDP),  and 


/?,(p)  =  w  RR»Ap)nXpA 7 

l  \  Ha  dt 


which  is  evaluated  along  the  Cagniard  contour.  p  is  the  ray  parameter,  17 a  is 
the  vertical  wave  slowness  (r]a  -  [a'2  -  p2]1/2  ),  RNZ(p)  is  the  receiver 
function  defined  by  Helmberger  (1974),  and  II1(p)  is  the  product  of  generalized 
reflection  and  transmission  coefficients  along  the  path  for  each  ray. 


The  form  of  the  RDP  is  generally  a  modification  of  that  proposed  by 
Haskell  (1967): 


r(()=  r.|  1  -p'*'  1  |  (3) 


where  'J',,  represents  the  static  value  with  dimensions  of  volume,  K  represents 
the  rise  time  or  corner  frequency  of  the  source  and  has  dimensions  of  inverse 


time,  and  B  represents  the  overshoot  or  the  extent  to  which  the  spectrum  is 
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peaked  at  its  corner,  and  is  dimensionless. 


von  Seggern  and  Blandford  (  1972.) 


truncated  the  polynomial  to  a  quadratic  based  on  the  observed  ^"2  spectra', 
decay  of  teleseismic,  short-period  body  waves.  While  Haskell's  source  model 
is  continuous  through  acceleration,  that  of  von  Seggern  and  Blandford  is 
continuous  onlv  to  displacement,  is  discontinuous  to  velocity  and  singular  to 
acceleration.  This  is  obviously  insufficient  for  modeling  local  velocities 
and  accelerations. 

As  a  compromise,  we  will  use  the  cubic  source  function  of  Helmberger  and 
Hadley  (1981): 


y(0-  1  -e~A< 


kt*^(kt)2- B(kt)3 


(•O 


which  is  continuous  through  velocity.  The  parameters  have  the  same  meaning 
between  source  models,  and  in  fact,  is  the  same  for  all  three  of  the  models 
discussed,  but  the  values  of  K  and  B  will  vary  between  source  models.  Barker, 
et  al.  (  1985b)  demonstrated  that  this  cubic  form  of  the  RDP ,  when  convolved 
with  elastic  Green's  functions,  can  explain  near-field  observations.  A 
quadratic  RDP  can  explain  the  observations  as  long  as  the  Green's  functions 
include  the  effects  of  depth-dependent  attenuation.  Rather  than  attempting  to 
solve  for  the  Q  structure  of  Pahute  Mesa,  we  have  chosen  the  Helmberger-Hadlev 
RDP  and  elastic  propagation  as  the  simpler  model.  One  disadvantage  is  that 
while  the  von  Seggern- Blandford  source  may  be  scaled  to  the  empirical  RDP  of 
Mueller  and  Murphy  (1971),  we  must  derive  scaling  relations  for  the 
Helm.berger-Hadley  source.  This  is  done  in  a  following  section.  As  discussed 
by  Burdick,  et  al .  (1982)  and  Hartzell,  et  al .  (1983),  the  short  durations  of 
locally  recorded  ground  motion  provide  poor  resolution  of  frequencies  well 
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below  the  source  corner,  so  B  trades  off  with  4'„ .  '’2  therefore  fix  B  -  1.0, 
and  solve  for  K  and  4>c .  Since  the  scaling  relations  are  derived  within  this 
constraint,  errors  in  the  assumption  of  B  will  not  adversely  affect  yield 
estimates,  but  should  be  taken  into  account  when  applied  to  lower  frequency 


The  results  of  forward  modeling  of  the  events  listed  in  Table  1  were 
presented  in  detail  in  Barker,  et  al .  (1985a)  and  are  summarized  in  Table  U. 
The  fits  of  the  synthetic  waveforms  to  the  observed  vertical  velocity  records 
from  outside  of  the  spall  zone  will  be  compared  with  the  results  of  the 
inversion  to  follow.  Barker,  et  al .  (1985a)  determined  a  scaling  relation 
using  the  model  parameters  and  K  to  estimate  the  yields  shown  in  Table  U , 
as  well  as  yields  for  four  additional  Pahute  Mesa  events. 


Table  U  -  Source  Parameters  from  Forward  Modeling 


Event 


BOXCAR 

MAST 

INLET 

SCOTCH 


(xlO10  cm3) 


Yield  published  Yield 

(kt)  (kt) 


1300  1300* 

943  5202 

286  3002 

127  155* 


:  Marshall,  et  al .  (1979) 

2  Dahlman  and  Israelson  (1977) 
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WAVEFORM  INVERSION  METHOD 


The  layered  structure  models  required  for  generalized  ray  synthetics  rrav 
be  considered  discretizations  of  one  or  more  velocity  gradients  separated  bv 
first-  or  second-order  discontinuities.  In  order  to  minimize  the  number  of 
free  parameters  in  the  inversion,  we  will  utilize  this  fact  and  solve  not  for 
the  velocities  and  thicknesses  of  individual  layers,  but  for  the  velocity 
gradient  and  the  velocity  at  the  top  of  the  gradient.  As  showm  in  Figure  6, 
the  discretized  model  follows  this  gradient,  with  layer  thicknesses  determined 
by  the  wavelength  of  the  predominant  period  in  the  data.  Shear  velocity  and 
density  may  be  free  parameters,  or  may  be  related  to  the  compressional 
velocity  by  predefined  constants  In  general,  the  free  parameters  are  the 
gradient  within  each  layer,  the  velocity  at  the  top  of  the  layer,  and  the 
depth  to  the  top  of  the  layer.  Thus,  for  an  n-layered  structure  in  which 
shear  velocity  and  density  are  not  free  parameters,  the  number  of  structure 
parameters  in  the  inversion  is  3n-l.  The  number  of  gradients  to  be  considered 
will  be  fixed  for  each  inversion,  so  the  optimal  number  will  have  to  be 
determined  by  test  inversions  proceeding  from  simpler  to  more  complex  models. 
Constraints  such  as  the  depth  to  the  water  table,  or  that  a  particular 
interface  may  have  only  a  second-order  discontinuity,  reduce  the  number  of 
free  parameters  while  increasing  non-linearity  of  the  inversion. 
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Since  the  discretized  velocity  model  will  vary  from  iteration  to 
iteration  as  well  as  with  the  model  perturbations  required  to  compute  the 
numerical  partials,  the  ray  set  from  which  the  synthetics  are  computed  will 
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Figure  6.  Structure  model  parameterization.  The  layered  structure  model  is 
discretized  from  one  or  more  velocity  gradients  defined  by  the  values  of  the 
gradient,  the  velocity  at  the  top  of  the  gradient  and  the  depth  to  the  top  of 
the  gradient  The  depth  increment  is  defined  as  the  product  of  the  P-wave 
velocity  and  the  predominant  period  of  the  data  (typically  0.2 b  sec).  S  wave 
velocity  and  density  are  linearly  related  to  the  P-wave  velocity  The  case 
shown  consists  of  two  gradients  separated  by  a  second-order  discontinuity. 
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also  vary.  Barker,  et  al.  (1985a)  found  chat  the  first  one  or  two  cycles  of 
vertical  velocity  observations  at  near-field  ranges  were  well  modeled  bv 
including  the  upgoing  direct  ray.  a  sum  of  diving  rays  that  constitute 
turning,  reflecting  and  critically  refracting  energy,  and  a  sum  of  rays  that 
depart  upward,  reflect  from  the  free  surface  and  follow  diving  ray  paths.  If 
radial  observations  are  to  be  modeled,  or  if  later  portions  of  the  waveform, 
are  considered,  P-to-S  conversions  at  the  free  surface  and  at  discontinuities 
should  be  included.  An  algorithm  to  generate  this  ray  set  for  each  new 
layered  structure  is  simple  to  implement,  and  computational  time  is  saved  by 
culling  rays  which  arrive  outside  the  inversion  time  window,  or  whose 
amplitude  will  be  below  some  predefined  cutoff  level. 

The  linearized  least- squares  inverse  may  be  obtained  by  solving 

Jc  '  =  A  '  J_£> '  ( 5 ) 

in  which  Ac'  is  the  residual  vector,  containing  the  differences  between 
observed  and  synthetic  data,  A'  is  the  matrix  of  partial  derivatives,  and  An 1 
is  the  desired  vector  of  parameter  changes.  In  this  method,  we  solve  for  the 
changes  in  the  model  parameters  from  those  of  the  starting  or  the  previous 
iteration  model.  The  resulting  error  analysis  gives  not  the  variances  of  the 
model  parameters,  but  the  variances  of  the  parameter  changes.  In  order  to 
determine  error  estimates  of  the  final  model  parameters  we  must  propagate 
these  variances  from  iteration  to  iteration  in  some  manner. 
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Shaw  and  Orcutt  (1985)  use  a  different  approach  to  solve  the  linearized 
inverse  (which  they  attribute  to  Parker,  unpublished  manuscript,  1981'. 
Rather  than  inverting  equation  (5),  they  add  A’p-'  (where  contains  re¬ 

starting  model)  to  both  sides  and  invert 

J  =  Jr  '  A  '  p_c, '  =  A  '  pp  *  Ajj  '  .  =  A  '  p  '  (  6  ) 

to  solve  directly  for  the  new  model  parameters,  p’.  The  benefit  of  (6),  which 
they  term  the  "jumping"  method,  is  that  rather  than  minimizing  the  parameter 
changes  and  obtaining  a  model  dependent  on  the  starting  model,  one  minimizes 
(by  some  appropriate  norm)  the  resulting  model  itself.  An  additional  benefit 
is  that  while  a  priori  constraints  on  the  model  parameters  (after  Jackson, 
1979)  may  be  desired,  it  is  often  difficult  to  transform  these  into 
constraints  on  the  parameter  changes.  While  we  desire  the  ability  to 
constrain  parameters  and  to  determine  error  estimates  of  the  parameters  rather 
than  the  parameter  changes,  we  have  a  certain  level  of  confidence  that  our 
starting  model  (obtained  from  forward  modeling)  is  reasonably  close  to  the 
correct  solution.  We  therefore  append  a  smoothing  constraint  onto  equatior, 
f6>  such  that 

Po'^P'  (  /  ) 

where-  A  is  a  diagonal  matrix  of  Lagrange  multipliers.  Thus  the  inversion 
attempts  to  solve  one  equation  for  a  new  set  of  parameters  while 
simultaneously  attempting  to  equate  those  parameters  to  the  starting  for 
previous  iteration)  model.  The  relative  weights  of  these  conflicting 
e.uu'ion.s  are  adjusted  via  the  Lagrange  multipliers.  In  practice,  we  use  a 
■ingle  Lagrange  multiplier,  so  A  is  a  constant. 
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Perhaps  the  simplest  determination  of  the  residual  vector  is  a 
po  int -by - po int  difference  of  observed  and  synthetic  seismograms  within  the 
time  window.  This  implies,  however,  that  neighboring  points  in  a  seismogram 
represent  independent  observations,  which  is  generally  not  the  case.  On  the 
other  hand,  since  each  seismogram  reflects  a  lagged  sum  of  rays  which  sample 
different  portions  of  the  structure,  the  resolution  of  parameters  varies 
through  the  time  window,  and  a  single  observation  for  each  seismogram, 
underestimates  the  information  contained  in  the  waveform.  Burdick  and  Me  liman 
(1976)  utilized  an  error,  or  residual,  function  based  on  the  normalized 
cross -correlation  coefficient  between  observed  and  synthetic  waveforms.  That 
is  , 


where 


e,  =  1  -  max  o,  ®  s , 


(8) 


6,(0 


o,(0 

y!  J  o?(f ) d( 


(9) 


and 


0(0 


s,(0 


(10) 


o,(t>  is  the  ith  observed  seismogram,  s^t)  is  the  ith  synthetic  seismogram  and 
a  denotes  c ross - cor re lat ion .  This  error  function  has  become  popular  because 
i  r.  is  sensitive  to  the  shapes  of  the  waveforms,  while  the  normalization  makes 
it  insensitive  to  absolute  amplitude,  and  measuring  the  maximum  of  the 
cross -correlation  makes  it  insensitive  to  absolute  time.  In  the  near- field 
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problem,  however,  absolute  amplitude  and  time  are  well  calibrated,  and  should 
be  fit  bv  the  resultant  model.  Therefore,  in  addition  to  (8),  our  residual 
will  include  for  each  seismogram  the  relative  difference  in  the  normalization 
factors , 


I  o;(t)dt~.  j  s*(t)dt 


(11) 


o?(t)dt 


which  is  a  measure  of  absolute  amplitude  residual,  and  the  time  lag  to  the 
maximum  of  the  cross  -  corre lat ion ,  t  ,  which  is  a  measure  of  absolute  time 
residual.  The  residual  vector  is  now  defined  as 


Jc'  =  .a  ,  ,e2.a2  ,t2 . en,am.t„. 


(12) 


for  m  seismograms.  The  objective  function  to  be  minimized  by  the  inverse  is 


r  =  p'ef-af-  r 


(13) 


and  the  RMS  measure  of  the  fit  of  the  synthetics  to  the  data  is  defined  as 
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For  the  inversion,  we  must  compute  the  partial  derivatives  of  the 
residual  functions  with  respect  to  the  individual  parameters  of  the  source  and 
velocity  structure  models.  For  a  particular  residual,  say  e, ,  and  parameter, 
sav  p  the  partial  derivative  is 


d  p  ,  o  e ,  d  u  ' 
d  p !  du. ' d  p  l 


(15) 


where  w'  is  the  velocity  response  of  equation  (2).  If  the  parameter  is  >I'I  or 
B,  the  RDP  is  linear  in  these  parameters  and  the  partial  is  easily  determined 
analytically.  However,  as  seen  in  equation  (4),  the  RDP  is  nonlinear  in  K 
Similarly,  the  partial  derivative  with  respect  to  a  structure  parameter,  say 
v.  ,  reduces  to  the  evaluation  of 


P)  =  dR,(p)  dp 

d  i  !  dp  di  ! 


(16) 


where  R,(t)  is  defined  in  equation  (1).  In  this  case,  the  partial  consists  of 
evaluating  the  partials  of  the  various  factors  in  Rx(p)  with  respect  to  the 
complex  ray  parameter  along  the  Cagniard  contour,  while  the  contour  itself 
varies  with  changes  in  the  velocity  parameter.  This  evaluation  must  be  made 
for  each  ray,  then  summed  to  obtain  the  partial  derivative. 


Previous  structure  inversions  based  on  waveform  modeling  have  invoked  a 
number  of  approximations  in  order  to  simplify  equation  (16).  For  example, 
Mellman  (1980)  utilized  a  modified  first-motion  approximation,  such  that 
rapidly  varying  terms  in  R^p)  are  evaluated  only  at  the  geometric  ray 
parameter.  Similarly,  Brown  (1982),  Shaw  (1983)  and  Chapman  and  Orcutt  (1985) 
invoke  the  WKBJ  approximation,  but  this  breaks  down  in  the  presence  of  large 
velocity  gradients  or  where  discontinuities  are  encountered  near  the  turning 
depth.  For  the  current  application,  the  logic  required  to  evaluate  the 
validity  of  these  approximations  for  each  ray  would  outweigh  the  effort  at 
computing  the  response  itself.  We  have  therefore  opted  to  calculate  the 
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partial  derivatives  numerically-  For  simpler  problems,  this  computati 


at  lor.  vo u.c 


c  e  v  t  a  i  n  i  v 


be  quite  inefficient.  However,  given  the  complexity  of  the 


near- field  wave  propagation  problem,  the  computational  cost  is  offset  bv  tr.o 
improved  accuracy  of  the  current  approach.  For  a  particular  residual,  say  e., 
and  parameter,  sav  p,,  the  numerical  partial  derivative  is  computed  by,  rather 
than  (131. 


do,  =  o,  p.  -  J p,  -  o,  pL 
•P,  -1 P, 


(17) 


where  £p,  is  the  perturbation  to  the  starting  model  parameter,  Pj  .  The 
partials  5a. /cp.  and  5tt/5p.  have  similar  form. 


Now  that  we  have  defined  the  form  of  the  residual  and  parameter  vectors 
and  the  partial  derivative  matrix,  we  may  redefine  the  matrices  in  (6)  as 


J  -  5  3J' 


e  -  !*: 


.4=5  MU  3 


(18) 


in  which  S  is  the  covariance  matrix  of  the  residuals  and  W  is  a  weighting 
matrix  of  the  parameter  changes.  For  simplicity,  we  assume  that  errors  are 
uncorrelated,  so  S  is  diagonal.  If  an  estimate  of  the  residual  variance  is 
known  (e.g. .  timing  errors),  that  value  is  used  in  S .  Otherwise,  the  residual 
variances  are  inversely  related  to  the  relative  confidence  in  the  residuals 
The  basic  purpose  of  the  diagonal  matrix  V  is  to  counter  the  dimensionality  of 
the  parameterization  of  source  and  structure.  Obviously  velocity  gradier.’ 
(with  values  on  the  order  of  1  km/s/km)  and  't*  (values  of  10;:  cm- )  should  i. 


be  weighted  equally.  Due  to  nonlinearity,  the  weighting  factors  are  also  used 
to  stabilize  the  parameter  changes .  For  example,  if  a  parameter  changes  to  an 
unreasonable  value  or  oscillates  about  some  value,  the  weight  of  that 
parameter  is  decreased  and  the  inversion  is  repeated.  This  introduces  a 
certain  degree  of  uncertainty  in  interpreting  the  variances  of  the  parameters, 
since  these  now  include  the  effects  of  the  somewhat  arbitrary  weights 
necessary  due  to  nonlinearity.  Although  convergence  criteria  may  be  utilized 
to  determine  when  to  stop  an  inversion,  we  have  not  implemented  such  at. 
approach,  but  simply  allow  additional  sets  of  five  iterations  until 
convergence  is  deemed  complete  or  unattainable .  Finally,  following  Jackson 
(1979),  a  priori  constraints  may  be  placed  on  the  parameter  changes  so  that, 
for  example,  the  velocity  cannot  be  negative,  or  so  that  source  parameters 
must  fall  within  a  range  of  reasonable  values. 

The  weighted  version  of  (6)  may  be  solved  by  a  linearized  generalized 
inverse  (Wiggins.  1972)  such  that 

p  =  r  j ' 1  r 7  j  (19) 

in  which  A  is  the  diagonal  matrix  of  non-zero  singular  values  of  A,  V  is  the 
corresponding  matrix  of  eigenvectors  which  spans  the  parameter  space,  and  U  is 
the  matrix  of  eigenvectors  that  spans  the  residual  space.  Since  the  problem 
is  non-linear,  the  inversion  is  iterative  and  minimizes  both  the  leas t - square s 
error,  |A  £  -  A|2,  and  the  length  of  the  parameter  vector,  |£l2.  Since  small 
singular  values  result  in  large  parameter  variances,  an  a  priori  cut-off  is 
specified,  below  which  the  singular  value  is  assumed  to  be  zero,  and  tin- 
dimension  of  the  parameter  spar--  is  truncated.  The  variances  of  t!;>- 


(yo) 


and  include  the  effects  of  the  parameter  weighting  matrix.  The  parameter 
trade-offs  that  result  from  singular  value  truncation  may  be  determined  by  at: 
inspection  of  the  parameter  resolution  matrix, 

«--rir  (21) 

which  reduces  to  the  identity  matrix  in  the  case  of  perfect  resolution.  If 
during  test  inversions,  the  parameter  resolution  is  deemed  unacceptable, 
adjustments  may  be  made  in  the  parameter  weights  or  the  singular  value 
cut-off,  and  the  inversion  may  be  repeated.  The  importance  of  particular 
observations  or  residuals  in  the  inversion  may  be  determined  from  the 
information  density  matrix, 


S  ~  CL7  (22) 

or  its  diagonal,  the  data  importance  vector.  This  may  indicate,  for  example, 
that  the  new  parameters  are  most  sensitive  to  the  absolute  amplitude  residual 
at  nearer  receivers,  and  to  the  absolute  time  residual  at  more  distant, 
receivers.  It  also  illustrates  the  importance  of  the  parameter  smoothing 
constraint  via  the  Lagrange  multipliers.  Thus  the  data  importance  vector  from 
test  inversions  may  suggest  modifications  to  the  observation  weights  or  to  the 
Lagrange  multiplier. 

One  distinct  disadvantage  of  our  implementation  of  the  jumping  method  for 
the  parameters  rather  than  the  usual  generalized  inverse  for  the  parameter 
changes  is.  in  the  treatment  of  parameters  when  eigenvalues  are  truncated.  I r, 
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the  generalized  inverse,  if  an  eigenvalue  is  truncated,  the  contribution  to 
the  parameter  changes  due  to  that  eigenvalue  is  zero.  Thus,  if  tht 
eigenvector  corresponding  to  a  truncated  eigenvalue  is  in  the  direction  of  a 
particular  parameter  change,  that  parameter  change  is  set  to  zero  and  the 
parameter  is  unchanged  in  that  iteration.  In  the  jumping  method,  this 
situation  would  result  in  the  new  parameter  value  being  set  to  zero.  This  is 
obviously  undesirable  for  parameters  such  as  <<e,  K  and  the  velocity  at  the  top 
of  a  layer . 

In  general,  the  new  parameter  value  at  a  given  iteration  is  a  linear 
combination  of  contributions  from  each  of  the  eigenvalues  not  truncated  during 
that  iteration.  The  extent  to  which  the  new  parameter  results  from  resolved 
eigenvalues  is  therefore  given  by  the  total  of  the  lengths  of  the  resolved 
eigenvectors  in  the  direction  of  that  parameter.  That  length  is  the  square 
root  of  the  diagonal  of  the  resolution  matrix  corresponding  to  the  parameter. 
To  avoid  the  undesirable  zeroing  of  parameters  described  above,  we  have 
implemented  a  smoothing  procedure  implemented  at  the  end  of  each  iteration. 
If  the  value  of  a  parameter  at  the  start  of  an  iteration  is  given  by  p  and 
the  new  parameter  value  at  the  end  of  that  iteration  is  p  ,  the  smoothed 
parameter  value  is  given  by 

p!s1  =  R„P,  +  \  1  "  R„Po,  (^) 

where  R1;  is  the  value  of  the  diagonal  of  the  parameter  resolution  matrix  for 
parameter  i.  This  is  obviously  a  nonlinear  constraint  imposed  upon  the 
i r. version,  and  has  the  effect  of  somewhat  duplicating  the  parameter  smoothing 


vith  the  Lagrange  multipliers.  As  such  it  slows  convergence  to  a  solutioj 
hut  has  the  desired  effect  that  the  parameters  of  that  solution  reta: 
reasonable  values. 
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TEST  OF  THE  INVERSION  METHOD 


r 


A  nuir.be r  of  tests  of  the  generalized  inverse  for  structure  parameters  or 
the  simultaneous  inverse  for  structure  and  source  parameters  were  presented 
ar.d  discussed  in  Barker,  et  al.  (  1985b,  1986).  In  particular,  those  reports 
illustrated  the  shapes  of  the  residual  functions  in  the  parameter  space  (i.e., 
the  locations  of  local  and  global  residual  minima),  and  discussed  the 
parameter  trade-offs  in  this  non-linear  problem.  Since  we  have  now 
implemented  the  jumping  method  of  inversion  (Shaw  and  Orcutt,  1985),  we  have 
repeated  one  of  the  test  inversions  of  Barker,  et  al  (  1986)  using  this 
technique . 


A  synthetic  data  set  was  generated  using  a  Helmberger-Hadley  (1981)  RDP 
and  a  two- layer  gradient  structure.  The  velocity  structure  was  based  on  the 
model  Hartzell,  et  al .  (1983)  obtained  for  Pahute  Mesa,  and  consists  of  a 
surface  velocity  of  2.5  km/s  and  gradient  of  0.83  km/s/km  in  the  top  layer. 
The  second  layer  is  separated  from  the  first  by  a  second-order  discontinuity 
at  a  depth  of  3.0  km  and  has  a  gradient  of  0.20  km/s/km  below  that.  For  the 
source  function,  we  set  -  1.0  x  10’*°  cm3,  k-7.0  s'1  and  B-=l .  Vertical 
component  velocity  synthetics  were  computed  for  a  source  at  1.0  km  depth  and 
for  ranges  from  2  to  12  km.  In  Test  «  4  of  Barker,  et  al.  (1986),  this  data 
set  was  first  inverted  for  the  parameters  of  the  source  and  of  the  upper 
layer.  Syntl  ;tic  data  at  greater  ranges  were  then  added  in  order  to  resolve 
the  lower  layer  parameters. 


For  the  current  test,  we  apply  the  jumping  method  inversion  to  tin 
original  data  set  and  solve  for  the  parameters  of  both  layers  as  well  as  th< 
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scarce.  There  are  six  parameters  in  the  inversion:  K,  the  upper  laver 

gradient  .  the  velocity  at  the  surface,  the  lov.-r  layer  gradient ,  and  the  dept:, 
to  the  top  of  the  lover  laver.  The  par.vrett  r  weighting  matrix  is  diar  I.'. 
1,1.  0.1,  0.1.  0.1.  1.0'  for  these  parameters,  respectively.  The  covariance 

matrix  of  the  observations  is  specified  ar:  diag(0.1,  1.0,  0.01)  for  tr.<- 

waveform,  amplitude  and  time  residuals,  respectively,  at  each  range.  The 
Lagrange  multiplier  for  parameter  smoothing  was  set  to  1.0,  and  the  eigenvalue 
truncation  parameter  was  chosen  to  allow  full  resolution.  Since  no 
eigenvalues  are  truncated,  the  solution  smoothing  applied  at  the  end  of  each 
iteration  has  no  effect  for  this  test.  An  arbitrary  starting  model  was 
specified  and  five  iterations  were  performed.  The  resulting  parameters  and 
their  standard  errors  (square  root  of  the  variance)  are  listed  in  Table  5  for 
each  iteration.  A  comparison  of  the  synthetic  waveforms  for  each  iteration 
with  the  synthetic  data  set  is  shown  in  Figure  7.  In  the  figure,  the  data  and 

synthetics  are  plotted  on  the  same  scales,  the  peak  absolute  amplitude  is 
written  to  the  right  of  each  trace,  and  the  inversion  time  window  is  indicated 
by  arrows.  The  data  importance  vectors,  including  the  relative  importance  of 
the  parameter  smoothing  equations,  are  plotted  in  Figure  8.  The  parameter 
resolution  matrix  for  each  iteration  is  simply  the  identity  matrix,  and  is  not 
plotted. 

It  is  clear  from  an  inspection  of  the  waveforms  (Figure  1)  that  the 
amplitude  and  timing  errors  of  the  starting  model  are  almost  entire  !v 
corrected  in  the  first  iteration.  The  following  iterations  simplv  fine  tur.e 
the  model  for  amplitudes  and  times,  while  attempting  to  fit  the  details  of  the 
v.iVe  forms  This  may  be  seen  in  the  data  importance  vectors  (Figure  8, 
Luring  the  first  iteration,  the  amplitude  and  time  residuals  are  tr<  : 
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Figure  ?.  Results  of  the  test  inversion  for  velocity  structure  Synthetic 
seismograms  are  shown  for  the  starting  model  and  five  iterations  at  e  a 
receiver  The  data  seismograms  (also  synthetic)  are  plotted  below  Ai  row-, 
above  the  data  waveforms  indicate  the  inversion  time  windows  and  numbers  to 
the  right  are  the  peak  trace  amplitudes. 
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Figure  8  The  data  Importance  vectors  for  the  test  inversion.  Shown  are  the 
diagonals  of  the  information  density  matrix  for  each  iteration  including  the 
relative  importance  of  the  parameter  smoothing  equations.  The  length  of  the 
ticV'.s  indicate  the  relative  importance  of  the  residuals  at  each  receiver  in 
the  determination  of  the  solution. 


important  in  constraining  the  solution  :i  the  seci-:/i  1  «  uiidi.  t :  • 
i-:  or*  ar.h  ot"  the  waveform  residuals  inert  ro-s  wh.ile  that  of  t  he  i  •  u  :• 

10  "idu. :1s  decreases  In  the  third  through  fifth  iterations,  the  mere  d  i  <•  :  u  • 
waveforms  are  important  as  deeper  structural  parameters  are  adjusted  Tr.- 
importance  of  the  time  residuals  remains  high  for  all  iteration-.,  cons  t  r  i  n  i :. : 
the  parameters  to  remain  near  the  correct  solution.  The  parameter  smoothing 
via  the  Lagrange  multiplier  has  its  largest  etfect  ir.  constraining  the  second 
parameter,  K.  This  is  because  although  we  have  specified  full  resolution,  the 
value  of  K  trades  off  with  4<c  and  may  varv  substantially.  The  parameter 
smoothing  prevents  radical  changes  in  K,  allowing  uriform  convergence 

The  convergence  of  the  inversion  parameters  listed  in  Table  5  to  the 
correct  model  parameters  (listed  at  the  bottom,  of  the  table)  indicates  the 
relative  resolution  of  the  parameters.  The  velocity  at  the  surface  has  nearlv 
reached  its  correct  value  after  only  one  iteration.  This  is  because,  for  the 
ranges  considered,  the  time  residual  is  most  sensitive  to  this  parameter.  The 
gradient  of  the  top  layer  converges  more  slowly.  4<„,  is  reduced  in  the  first 
iteration  to  smaller  than  the  correct  value,  rebounds  in  the  second  iteration, 
the:,  converges  to  a  value  somewhat  below  the  correct  value.  K  decreases 
men."  topically  due  to  the  parameter  smoothing  constraint,  but  because  of  its 
t  rade-of i  with  .  converges  to  a  value  somewhat  too  high  For  these  ranges . 

he  r-sidu.iis  are  not  very  sensitive  to  the  gradient  in  the  lower  lave  r .  Ir 
:  u't,  it  appears  that  the  lower  gradient  and  d< pth  parameters  are  given  value- 
•  hat  compensate  for  the  upper  gradient  being  slightly  too  low.  As  in  Far her 
ei  a  1  ' 1 h  8  h  1  .  the  inclusion  of  data  from  greater  ranges  would  improve  th- 

resolution  of  the  lower  laver  parameters  Given  the  inherent  trade-off; 
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We  begin  the  inversions  of  observed  data  with  the  Par  ute  Mesa  ever. r. 
E LX CAR  v-  '26/  68.  1300  kt  announced,  1.160  km  depth).  Hartzell  et  al  ,  (19-;. 
developed  their  Pahute  Mesa  structure  model  primarily  by  modeling  the 
waveforms  of  receivers  S-24  and  S-34  from  BOXCAR,  although  travel  times  from, 
several  other  events  were  quite  well  predicted.  We  will  use.  an  approximation 
of  the  Hartzell,  et  al.  (1983)  model  as  our  starting  model.  In  addition, 
however,  we  wish  to  include  a  high  velocity  gradient  above  the  water  table  as 
suggested  by  the  velocity  log  published  in  Stump  and  Johnson  (1984).  This  drv 
tuff  layer  is  adequately  approximated  by  a  gradient  of  0.92  km/s/km  starting 
from  a  surface  velocity  of  2.4  km/s  and  reaching  the  water  table  at  a  depth  of 
0.65  km.  These  parameters  are  held  fixed  through  the  inversion.  In  the 
inversion,  we  allow  only  second-order  discontinuities  below  the  water  table, 
so  the  free  parameters  consist  of  K,  the  gradient  in  the  second  layer,  the 
gradient  in  the  third  layer,  and  the  depth  to  the  top  of  the  third  layer.  The 
parameter  weighting  matrix  is  specified  as  diag(10.,  10.,  1.0,  0.1,  0.1)  for 
these  parameters,  respectively.  The  starting  model  source  parameters  are  the 
forward  modeling  results  listed  in  Table  4. 

The  stations  used  in  the  inversion  are  listed  in  Table  2,  and  include  a 
broad  distribution  of  ranges.  The  closest  receivers  ( S  - 1 2  at  3.81  km  and  S - 1 6 
at  4.8?  km)  are  only  marginally  outside  the  spall  zone  and  appear  to  include 
arrivals  attributable  to  spall  opening  or  slap-down.  Since  these  effects  will 
not  be  included  in  our  modeling,  we  must  increase  the  covariance  of  the 
waveform  residuals,  effectively  down-weighting  the  importance  of  fitting  the 
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wave- form .  at  these  receivers 


The  covariance  matrix  is  diag(10.0,  10.0,  0.01) 


'  V 


f. 


V, 


lor  waveform,  amplitude  and  travel  time,  respectively.  Similarly,  the 
farthest  station  ( S - 7  a  at  22. km)  is  located  outside  the  Silent  Canvon 
Caldera  of  Pahute  Mesa,  and  we  could  expect  the  plane  -  layered  structure 
approximation  to  break  down  for  this  receiver.  We  have  therefore  increased  the 
covariance  of  the  waveform  and  travel-time  residuals  at  this  receiver.  The 
resulting  covariance  is  diag(10.0,  10.0,  10.0).  This  leaves  two  receivers 
(S-24  at  7.27  km  and  S-34  at  10.37  km)  with  reasonable  covariance  in  waveform, 
and  amplitude  residuals.  Unfortunately,  preliminary  inversions  indicate 
that  these  stations  have  unmodeled  travel-time  residuals  which  may  result  from, 
topography  or  laterally  varying  near-surface  velocities.  Accounting  for 


> 

j’- 


timing  errors, 
Five  iterations 
and  the  eigenva 


the  covariance  for  these  receivers  is  diag(0.1,  0.1 

were  allowed,  with  the  Lagrange  multiplier  specified 
lue  cutoff  set  for  full  resolution. 


1.0). 

as  1.0 


►  < 


The  inversion  results  are  listed  in  Table  6  and  the  waveforms  for  the 
final  inversion  model  are  plotted  in  Figure  9.  Also  shown  are  waveforms  that 
resulted  from  the  forward  modeling  study  of  Barker,  et  al .  (1985a)  using  the 
structure  model  in  Table  3  and  the  source  parameters  in  Table  U.  For  each 
station,  the  waveforms  are  plotted  on  the  same  time  and  amplitude  scales,  and 
the  peak  absolute  amplitude  is  written  to  the  right  of  each  trace.  Arrows 
above  the  observed  records  indicate  the  inversion  time  window  used.  In  both 
the  forward  and  inverse  synthetics,  the  waveforms  at  S - 1 2  and  S-16  are  simple 
and  impulsive,  while  the  observed  waveforms  are  broad,  and  include  significant 
late  arrivals.  These  may  be  attributed  to  unmodeled  processes  of  spall  and 
slap-down.  The  first  peak  at  S-24  and  S-34  includes  the  constructive 
interference  of  upgoirg  and  diving  rays,  while  the  trough  includes  the 
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able  6  -  BOXCAR  Inversion  Results 


e rat  i  on 

(,xl0lc  cm3) 

k 

(s'1) 

aaz/dz 
(km/s  /km. ) 

ca  ,  /  dz 
( km/s /km ) 

Z3 

(km) 

RMS 

fit 

LSE 

Start 

12.0 

6.50 

0.83 

0.20 

3.00 

.247 

l 

12  .  3 

6  .  56 

0.80 

0.21 

2 . 99 

.299 

0.90 

(1.7) 

(1.16) 

(.26) 

(.36) 

(  .35) 

2 

12 . 6 

6.62 

0.83 

0.18 

2 . 99 

.253 

0.98 

(1.7) 

(1.13) 

(29) 

(-29) 

(  .56) 

3 

12.8 

6.62 

0.82 

0.18 

2 . 99 

.  286 

0.91 

(1.7) 
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(  .34) 

(.37) 

(.56) 
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C  .  15 
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(1.18) 

(.27) 

(  ■  30) 

(.55) 

5 

13  .  2 

6.53 

0.81 

0.14 

2.99 

.237 

1.00 

(1.7) 

(1.06) 

(  .25) 

(.31) 

(.56) 

LSE  is  the  leas t - squares  error. 

Numbers  in  parentheses  are  the  standard  errors  of  the  parameters. 

Fixed  structure  parameters  include  a,_-2.4  km/s ,  Sc^-0.92  km/s/km  and 

zz-0 .65  km . 


Modeling  Results  for  BOXCAR 
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Figure  9.  Inversion  results  for  Pahute  Mesa  event  BOXCAR.  Shown  are  the 
observed  records  (top),  the  results  from  the  numerical  inversion  (center)  and 
from  forward  modeling  (bottom).  The  peak  trace  amplitude  is  shown  to  the 
right  and  the  inversion  time  window  is  denoted  by  arrows. 
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imre” action  of  the  surface  reflection  phases  Although  the  synthetic  arrival 
riit.es  at  these  stations  is  somewhat  off  iS-29  is  late  bv  0.02  s  while  S  -  2  *  is 
early  hy  0.09  si .  the  details  of  the  waveforms  are  quite  well  modeled  both  ir. 
the  inversion,  and  in  the  forward  modeling,  indicating  that  the  velocity 
gradient  near  the  source  depth  is  well  determined.  Since  S  -  79  is  outside  of 
Tribute  Mesa,  we  expect  the  arrival  time  to  be  poorlv  modeled.  The  amplitudes 
of  the  synthetics  are  also  somewhat  too  large,  but  the  general  shape  of  the 
waveforms  is  adequately  modeled  within  the  inversion  time  window. 

Table  6  indicates  that  in  the  inversion,  the  model  parameters  vary  or.lv 
slightly  from  the  starting  model.  K  increases,  but  then  returns  to  its 
starting  value.  <<„.  increases  slowly.  The  second  layer  gradient  varies  about 
its  starting  value,  while  the  third  layer  gradient  decreases  slowly  and  the 
depth  to  the  third  layer  remains  constant.  This  is  reflected  in  the  data 
importance  vectors  (Figure  10).  Throughout  the  inversion  the  amplitude 
residuals  (particularly  at  the  more  distant  stations)  are  most  important  it. 
determining  the  solution,  while  the  parameter  smoothing  constraint  on  iy  and 
the  third  layer  depth  force  these  parameters  to  vary  slowly.  Although  full 
resolution  was  specified,  the  problem  is  non-linear  and  the  parameters  trade 
off.  The  parameter  smoothing  inhibits  the  parameters  from  assuming 
unrealistic  values.  A  side  effect  of  full  resolution  is  that  the  parameter 
variances  can  be  large.  Shown  in  parentheses  in  Table  6  are  the  standard 
errors  (square  root  of  the  variances)  for  the  parameters.  Although  the  errors 
'-‘•'e  large  (particularly  for  the  third  layer  gradient),  the  solution  is 
c  or w  rgent  . 
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Smoothing  Time  Amplitude  Waveform 
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sr.ou .  a .  pernaps ,  ne  no  surprise  mat  me  inversion  re-sv-ts  co  n 
vie:  ir:  significantly  from  the  starting  model  .  The  starting  structure  model  i  5. 
:  used  or.  the  velocity  structure  determined  tv  Hartzell ,  et  al  .  (  1983  )  larg. 

:  rom  n  .•'deling  the  details  of  the  S-  24  and  S  -  3-  records  from  BOXCAR .  This  :.a< 
hi-::  si: ovr.  in  Figures  3  and  4,  and  is  reinforced  bv  the  agreement  between 
forward  and  inverse  synthetics  in  Figure  9.  The  conclusion  is  that  this, 
structure  model  is  very  good  for  BOXCAR.  The  starting  model  source  parameters 
were  determined  by  forward  modeling  using  this  structure  model.  One  problem, 
with  the  forward  modeling  study  was  the  need  to  assume  the  same  structure 
model  for  all  Pahute  Mesa  events.  By  performing  the  simultaneous  inversion  on 
the  waveforms  of  these  events,  we  will  determine  if  significant  modi f icat ions 
to  the  structure  model  are  demanded  for  slightly  different  locations  within 
Pahute  Mesa.  By  comparing  the  source  parameters  from  forward  and  inverse 
modeling,  we  may  determine  to  what  extent  these  structure  variations  effect 
our  estimates  of  source  parameters. 
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MA  S  T 

In  inverting  the  vertical  velocity  observations  for  the  event  MAST 
(6/19/75,  520  kt  estimated,  0.912  km  depth),  we  use  a  starting  model  with  the 
same  structure  parameters  as  in  the  BOXCAR  inversion  and  source  parameters 
determined  bv  forward  modeling  of  MAST  (Table  9).  Three  observed  records  are 
available  outside  the  spall  zone  (Table  2).  The  covariance  matrix  of  the 
observations  is  specified  as  diag(0.01,  0.1,  0.01)  for  waveform,  amplitude  and 
time  residuals,  respectively,  except  that  because  of  noise  at  the  most  distant 
receiver,  its  waveform  residual  variance  is  set  to  0.1.  The  parameter 
weighting  matrix  is  diag(1.0,  1.0,  0.1,  0.1,  0.1)  for  K,  the  gradients  in 
the  second  and  third  layers,  and  the  depth  to  the  top  of  the  third  layer, 
respectively.  The  Lagrange  multiplier  is  1.0  and  the  resolution  is  such  that 
three  or  four  of  the  five  parameters  are  resolved  each  iteration. 


The  results  of  the  inversion  are  listed  in  Table  7.  The  synthetic 
waveforms  from  the  fourth  iteration  are  plotted  in  Figure  11,  along  with  the 
observed  records  and  the  synthetics  from  the  forward  modeling  study  of  Barker, 
et  al.  '1985a).  Arrows  above  the  observed  data  record  indicate  the  inversion 
time  window,  and  the  peak  absolute  amplitude  of  each  trace  is  written  to  the 
right  The  parameter  resolution  matrices  for  each  iteration  are  plotted  in 
Figure  12.  Because  some  eigenvaluts  are  truncated,  the  gradient  in  the  third 
layer  is  not  resolved  and  other  parameters  trade  off  with  each  other  (in  the 
case  of  perfect  resolution,  the  parameter  resolution  matrix  would  be  the 
identity  matrix  i  K  and  the  second  layer  gradient  are  well  resolved.  is 
initially  well  resolved,  but  begins  to  trade  off  with  the  third  layer 
parameters  beginning  with  the  third  iteration.  The  inversion  reaches  a 
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LSE  is  the  least  -  squares  error. 

Numbers  in  parentheses  are  the  standard  errors  of  the  parameters. 

Fixed  structure  parameters  include  a. -2.  U  km/s ,  cq.-0.92  km/s/km  ar.c 

c,=0  65  km. 


Modeling  Results  for  MAST 


Observed  Data 
Inversion  Result 
Forward  Result 


S-5  3.65  km 


43.51  cm/s 
55.22  cm/s 
45.37  cm/s 


Observed  Data 
Inversion  Result 
Forward  Result 

Observed  Data 
Inversion  Result 
Forward  Result 


S-6  5.47  km 


yv 


S-7  7.30  km 


Time  (sec) 


19.36  cm/s 
18.08  cm/s 
33.24  cm/s 


13.73  cm/s 
12.06  cm/s 
17.72  cm/s 


Figure  11  Inversion  results  for  Pahute  Mesa  event  MAST.  Shown  are 
observed  records  (top),  the  results  from  the  numerical  inversion  (center) 
from  forward  modeling  (bottom)  The  peak,  trace  amplitude  is  shown  to 
right  and  the  inversion  time  window  is  denoted  by  arrows. 
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Figure  12.  Parameter  resolution  matrices  for  each  iteratio 
inversion  Plotted  from  top  to  bottom  are  the  rows  of 
corresponding  to  the  source  parameters,  ».  and  k.  and 
parameters ,  upper  gradient,  lover  gradient  and  interface  dept! 
resolution,  the  resolution  matrix  would  be  tht  identity  matrix 
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minimum  in  RMS  fit  at  the  fourth  iteration  (the  divergence  at  the  fifth 
iteration  is  due  primarily  to  the  change  in  'I'*),  so  the  synthetics  plotted  in. 
figure  11  represent  this  solution.  Presumably,  if  more  iterations  wv  rv 
allowed,  the  inversion  would  return  toward  these  parameters.  Because 
eigenvalue  truncation  was  allowed.  the  standard  errors  in  Table  t  are 
reasonably  small.  The  parameter  variance  has  improved  at  the  expense  of 
raram.eter  resolution. 


The  first  peak  at  S-5  is  simple  and  is  well  modeled  bv  both  the  forward 
ar.d  inverse  models.  The  trough  of  this  record.  howv'.vr,  reflects  the 
interference  of  the  pP  phase  with  the  upgoing  and  diving  P  waves.  The  details 
of  this  interference  are  slightly  better  fit  by  the  inversion  waveform,  but 
the  first  peak  amplitude  is  better  matched  by  the  forward  model.  The  features 
of  the  first  peak  at  S-6  and  S-7  are  sensitive  to  the  interference  of  upgoing 
ar.d  diving  P , and  are  significantly  better  modeled  by  the  inversion  waveforms. 
Further,  the  details  within  the  troughs  for  these  stations,  the  amplitudes  of 
the  first  peaks,  and  the  travel  times  are  all  improved  bv  the  inversion.  It 
would  appear  that  the  higher  velocity  gradient  in  the  second  layer  resolved  by 
the  inversion  is  demanded  by  the  data.  This  results  in  a  source  model  in 
which  is  somewhat  higher  and  K  is  somewhat  lower  than  the  forward  modeling 
results  of  Barker,  et  al.  (1985a). 
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INLET 


The  same  starting  structure  model  is  usee  to  invert  the  observed  vertical 
velocity  waveforms  for  the  INLET  (11/20/75,  500  Kt  estimated,  0.817km  depth) , 
and  starting  source  parameters  appropriate  for  this  event  are  obtained  fret 
the  forward  modeling  results  (Table  4).  Three  observed  records  are  available 
outside  the  spall  zone  (Table  2).  Although  the  spall  zone  for  INLET  was  quite- 
limited,  the  record  at  1.63  km  may  contain  some  unmodeled  effects  due  to  spall 
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Livers,  and  the  depth  to  the  top  of  the  third  layer,  respectively.  The 
Lagrange  multiplier  is  1.0  and  the  resolution  is  such  that  three  or  four  of 
the  five  parameters  are  resolved  each  iteration. 

The  results  of  the  inversion  are  listed  in  Table  8  The  synthetic 
waveforms  are  plotted  in  Figure  13,  along  with  the  observed  records  and  the 
r  vuthet  ics  from  the  forward  modeling  study  of  Barker,  et  al  ,  (1985a-' .  Arrows 

.■hove  the  observed  data  indicate  the  inversion  time  window,  and  the  peal, 
absolute  amplitude  of  each  trace  is  written,  to  the  right.  The  synthetic 
w.iVe  f  orn  from  forward  and  inverse  modeling  are  quite  similar  in  shape. 
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Table  8 


INLET  Inversion  Results 


I  teration 


(xlO--  crrr.) 


k 

(s’1) 


Za2,/Zz 
(kit  •  s  /km) 


Za3  ,/Zz 
(km/s  /km1) 


( km ") 


RMS 

fit 


LSE 


Modeling  Results  for  INLET 


Observed  Data 
Inversion  Result 
Forward  Result 


S-5  1.63  km 


95.40  cm/s 
83.69  cm/s 
78.38  cm/s 


Observed  Data 
Inversion  Result 
Forward  Result 


Observed  Data 
Inversion  Result 
Forward  Result 


S-6  3.27  km 


S-7  6.53  km 


40.82  cm/s 
40.97  cm/s 
39.29  cm/s 

22.07  cm/s 
24.57  cm/s 
19.90  cm/s 


0.0  1.0 

Time  (sec) 


~ i 
2.0 


Figure  13.  Inversion  results  for  Pahute  Mesa  event  INLET  Shown  are  t 
observed  records  (top),  the  results  from  the  numerical  inversion  (center)  a 
from  forward  modeling  (bottom).  The  peak  trace  amplitude  is  shown,  to  t 
right  and  the  inversion  time  window  is  denoted  by  arrows 
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although  the  absolute  amplitudes  and  travel  times  are  slightly  improved  bv  the 


inversion.  The  S-5  observed  record  contains  late  arrivals  in  the  trough  that 
are  not  modeled  by  either  approach  and  may  reflect  spall  processes.  The 
observed  first  pulse  at  S-6  is  significantly  wider  than  either  of  the 
synthetics,  although  the  pulsevidth  at  S-7  is  well  modeled  and  that  at  S-5  is 
nearly  fit.  Apparently,  allowing  the  structure  model  to  varv  cannot  improve 
the  fit  at  S-6  while  still  fitting  S-5  and  S-7.  The  observed  waveform  at  S-~ 
is  rather  noisy,  but  both  the  forward  and  inverse  synthetics  include  the 
multiple  arrivals  apparent  in  the  trough. 


The  parameter  resolution  matrix  for  INLET  is  plotted  in  Figure  14.  The 
gradient  in  the  third  layer  is  never  resolved.  The  slight  variations  in  this 
parameter  from  the  starting  value  (Table  8)  are  not  resolved,  but  result  from 
slight  trade-offs  with  other  parameters.  The  gradient  in  the  second  layer  is 
well  resolved  in  each  iteration,  but  its  value  departs  only  slightly  from  that 
of  the  starting  model.  This  indicates  that  the  velocity  gradient  near  the 
source  depth  from  the  Hartzell,  et  al .  (1983)  model  is  appropriate  for  INLET 

as  well  as  for  BOXCAR.  In  the  inversion,  K  is  well  resolved,  but  remains 
close  to  the  forward  modeling  value  determined  by  Barker,  et  al.  (1985a). 
however,  increases  by  about  35%  over  the  forward  modeling  value.  This  results 
in  the  improved  fits  to  the  peak  amplitudes. 


I  ^  < 


SCOTCH 


For  the  Pahute  Mesa  event  S  1 CH  ( 5/2  3  .'6  7  .  155  Kt  announced,  0 . 9’0  kr 
depth),  onlv  two  vertical  velocity  recordings  are  available  outside  the  spall 
.tone  but  within  Pahute  Mesa  (  see  Table  2)  Since  the  core  distant  of  these  is 
at  a  range  of  onlv  6.06  km,  and  i  rorr.  our  experience  in  inverting  MAST  and 
INLET  data,  we  suspect  that  the  parameters  of  the  lower  structure  will  be 
pocrlv  resolved.  We  therefore  use  the  same  starting  model  as  in  the  other 
events,  but  raise  the  interface  at  the  top  of  the  third  layer  from  3.0  km  to 
13  km.  In  addition,  we  decrease  the  weights  of  these  parameters.  The 
starting  source  model  is  obtained  from  the  forward  modeling  results  (Table  4). 
The  parameter  weighting  matrix  is  diag(1.0,  1.0,  0.1,  0.01,  0.01)  for  4'.,  K, 
the  gradients  of  the  second  and  third  layers,  and  the  depth  to  the  third 
layer,  respectively.  The  covariance  matrix  of  the  observations  is  diag(0.01, 
0.1,  0.001)  for  waveform,  amplitude  and  time  residuals,  respectively,  for  both 
stations.  The  Lagrange  multiplier  is,  once  again,  1.0,  and  the  resolution  is 
such  that  four  of  the  five  eigenvalues  are  resolved. 
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The  results  of  the  inversion  are  listed  in  Table  9,  and  the  final 
iteration  synthetic  seismograms  are  plotted  in  Figure  15  along  with  the 
observed  waveforms  and  the  synthetics  from  forward  modeling.  The  parameter 
resolution  matrices  for  each  iteration  are  plotted  in  Figure  16.  As  expected, 
the  gradient  in  rite  third  layer  is  not  resolved  in  the  early  iterations,  but 
because  the  eigenvector  corresponding  to  the  truncated  eigenvalue  is  not 
entirely  in  the  direction  of  this  parameter,  its  value  is  decreased  in  spite 
of  tie-  parameter  smoothing.  An  interesting  effect  occurs  at  the  fourth 
iteration,  however  At  this  point,  the  depth  to  the  top  of  the  third  layer 


Table  9  -  SCOTCH  Inversion  Results 


reratio 

n 


( x 1 0 - '  cm3)  (s'1) 


da2/d z  3a y  ?z  z,  RMS  LSE 

(km/s/km  (km /s  km  (km)  fit 

)  I 


3  .  35 

11 . 55 

1  .  10 

0 .63 

0 .92 

.451  162. 

(  .44) 

(  .87) 

(.20) 

(  .28) 

(  -15) 

2.02 

10.50 

1.00 

0.63 

0.76 

.229  20.0 

(  .62) 

(.88) 

(.12) 

(  .24) 

(  .28) 

LSE  is  the  leas t - squares  error. 

Numbers  in  parentheses  are  the  standard  errors  of  the  parameters. 

Fixed  structure  parameters  include  a-2  . 4  k m/s,  oa--0. 92  km/s/km  and 
.65  km . 


Modeling  Results  for  SCOTCH 


Figure  15.  Inversion  results  for  Pahute  Mesa  event  SCOTCH.  Shown  are  the 
observed  records  (top),  the  results  from  the  numerical  inversion  (center)  and 
frorr.  forward  modeling  (bottom).  The  peak  trace  amplitude  is  shown  to  the 
right  and  the  inversion  time  window  is  denoted  by  arrows. 
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SCOTCH  Resolution  Matrices 


Iteration 

Figure  16,  Parameter  resolution  matrices  for  each  iteration  of  the  SCOTCH 
inversion.  Plotted  from  top  to  bottom  are  the  rows  of  the  matrices 
corresponding  to  the  source  parameters,  f.  and  k,  and  the  structure 
parameters,  upper  gradient,  lower  gradient  and  interface  depth.  For  perfect 
resolution,  the  resolution  matrix  would  be  the  identity  matrix. 
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brcorr.es  shallower  (0  915  km)  than  the  source  depth  (0.970  km.)  ,  so  the  source 
it.  located  within  the  third  layer .  Now  the  gradient  in  the  third  laver  is 
well  resolved,  and  approaches  the  value  previously  obtained  for  the  second 
laver  (about  0.6  k:r.  s  'km),  while  the  second  layer  gradient  approaches  t hr 
value  of  the  fined  gradient  in  the  top  layer  (about  0.9  km/ s /kail  .  Sv  the 
fifth  iteration,  the  depth  to  the  top  of  the  third  layer  has  decreased  to  0." 
k.T.  The  result  of  all  of  these  changes  is  that  the  velocity  structure 
approximates  a  two-gradient  model  with  a  second  order  di scont inui t v  at  a  depth, 
of  0.77  km,  rather  than  at  the  presumed  water  table  depth  of  0.65  km  (we  do 
not  have  specific  information  on  the  water  table  depth  near  SCOTCH). 

The  observed  waveforms  (Figure  15)  include  an  impulsive  first  arrival.  A 


second  phase  arrives  at  a  time  between  the  first  peak  and  the  trough.  These 
features  are  well  modeled  by  the  forward  synthetic  seismogram  at  S-3A.  but  the 
two  arrivals  form  a  single  peak  in.  the  synthetic  for  S-9,  which  is  broader 
than  the  observed  pulse  Because  the  gradient  continues  below  the  source  in 
the  inversion  model,  the  diving  ray  arrives  before  the  direct  ray  for  both 
ranges  in  the  inversion  synthetics.  This  causes  a  somewhat  early  and 
doub 1 e - peaked  first  pulse.  Forcing  the  structure  below  the  source  to  have  a 
lower  gradient  would  improve  the  fit  of  the  synthetics  to  the  data. 
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SCALING  RELATION  FOR  YIELD  ESTIMATION 


The  source  parameters  determined  through  waveform  inversion  for  the  four 
Fahute  Mesa  events  studied  are  summarized  in  Table  10  Since  we  have  chosen 
an  analytic  form  of  the  RDP  (the  Helmberger  -  Hadlev ,  1981,  form),  we  must 
derive  scaling  relations  in  order  to  estimate  the  Yields  of  these  events. 
Burdick,  et  al  .  (  1982)  and  Hartzell,  et  al .  (  1983)  determined  scaling 
relations  for  the  Helmberger -Hadley  RDP  appropriate  for  Arr.chitka  and  Pahute 
Mesa,  respectively.  However,  each  derived  separate  relations  for  and  K  as 
functions  of  yield  and  source  depth.  Since  the  determinations  of  <l'c  ar.d  K  are 
not  independent,  we  follow  Barker,  et  al .  (1985a)  and  equate  the  relations  bv 
solving  for  the  common  parameter  of  depth,  to  obtain  a  single  scaling  relation 
giving  yield  as  a  function  of  and  K.  Specifically,  the  general  form  of  the 
relation  between  yield,  'I'*,  and  depth  is 


>  =C,Vl'h 


(?4) 


where  Cj  is  a  proportionality  constant,  and  a  and  b  are  exponents  of  yield  and 
depth,  respectively,  by  which  cavity  radius  is  related.  Similarly,  the 
relation  between  yield,  K  and  depth  has  the  form 


>'  =  C- 


(2‘o) 


where  C2  is  another  proportionality  constant,  and  n  is  an  empirically  derived 
constant.  Solving  these  equations  for  depth  and  equating  the  two,  we  obtain 


)  =  3mm  ^  ’• 


(26) 


Murphy  and  Mueller  (1971)  argue  that  it  is  reasonable  for  the  cavity  radius  to 
depend  on  the  cube  root  of  yield,  so  a  -  1/3.  Using  Orphal's  (1970)  equation 
for  cavity  radius,  we  obtain  b  -  0.09.  Finally,  Mueller  and  Murphy  ( 1 9  ~  1  "■ 
determined  empirically  that  n  -  2.4,  so  that  equation  (26)  becomes 

>  =  r,  v'l 7?b  k G  826  (?/) 

Of  the  four  events  modeled,  BOXCAR  and  SCOTCH  have  announced  yields,  and 
so  mav  be  used  to  determine  the  proportionality  constant,  C-, .  Using  the 
parameters  of  the  RDP  determined  for  BOXCAR  (Table  10)  and  a  vield  of  1300  kt 
(Marshall,  et  al . ,  1979),  we  obtain  C3  -  1.779  x  10"12.  This  value  predicts  a 
yield  of  175  kt  for  SCOTCH,  which  is  greater  than  the  announced  vield  of  155 
kt  (Marshall,  et  al.,  1979).  If,  on  the  other  hand,  we  use  the  yield  and  RDP 
parameters  for  SCOTCH  (Table  10),  ve  obtain  C3  -  1.575  x  10'12  which 
underpredicts  the  yield  of  BOXCAR  (1150  kt)  .  Predicted  yields  using  both 
proportionality  constants  are  given  in  Table  10,  but  since  SCOTCH  is  a  small, 
overburied  event,  we  prefer  the  predictions  based  on  the  BOXC.AR  constant. 


Table  10  -  Source  Parameters  from  Wave  form  Inversion 


Event 

K 

( s  ~ : ) 

4i 

*  cc 

(x!0ic  cm3) 

Yield- 

(kt) 

Yie Id2 
(kt) 

publ i shed 
(kt) 

BOXCAR 

6  .  53 

13.22 

1300* 

1150 

13003 

MAST 

6  .  39 

6.22 

497 

432 

520* 

inlet 

8  .  12 

4.48 

391 

346 

500* 

SCOTCH 

10. 50 

2.02 

175 

155* 

1 5  5 3 

Constrained  to  announced  yield. 

•  From  equation  (27)  calibrated  to  BOXCAR  ( C ,  -  1.779  x  10"*2) 

2  From  equation  (27)  calibrated  to  SCOTCH  (C,  -  1.575  x  10~uj 

3  Marshall,  et  al.  (  1 9 7 9 ) 

*  Dahl  man  and  Israelson  (1977) 
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CONCLUSIONS 


Perhaps  the  most  crucial  step  in  the  determination  of  effective  source 
functions  and  the  estimation  of  yields  from  near-field  recordings  of 
underground  nuclear  explosions  is  the  development  of  an  accurate  velocity- 
structure  "odel.  Previous  studies  have  obtained  these  models  through 
laborious  trial - and-error  waveform  modeling.  We  have  reviewed  the  results  of 
forward  modeling  for  source  parameters  at  Pahute  Mesa  using  the  velocity- 
structure  model  of  Hartzell,  et  al  .  (  1983).  The  resulting  synthetic 
seismograms  provide  an  admirable  fit  to  the  observed  waveforms,  particularly 
f->r  PCXCAP,  and  SCOTCH.  There  is  some  concern,  however,  regarding  the  adequacy 
■of  the  structure  model  for  other  locations  within  Pahute  Mesa,  and  in  how 
errors  in  the  structure  model  effect  the  determination  of  source  parameters. 
An  additional  problem  is  that  in  determining  source  parameters  from  explosions 
in  a  different  structure  (whether  Yucca  Flats,  Amchitka,  or  one  of  the  Soviet 
test  sites),  a  new  velocity  structure  model  must  first  be  determined. 

In  order  to  address  these  issues,  we  have  developed  a  simultaneous 
inversion  method  which  solves  for  the  parameters  of  both  the  source  and 
velocity  structure  models.  A  report  on  the  initial  development  of  the  method 
was  given  in  Barker,  et  al.  (198th)  and  the  application  of  this  method  to 
Pahute  Mesa  data  was  presented  in  Barker,  et  al .  (1986).  We  have  improved  the 
inversion  technique  by  incorporating  the  jumping  method  of  Shaw  and  Orcutt 
(1983).  The  primary  improvement  in  this  technique  is  that  parameter 
constraints  and  error  analyses  are  made  on  the  parameters  of  the  model,  rather 
than  on  the  changes  in  the  parameters.  We  have  presented  a  test  inversion  of 


a  svnthetic  data  set.  The  results  were  quit-.-  promising,  but  also  illustrate  : 
the  1  irritations  in  the  resolution  of  deeper  structures  when  limited  ranges  :  • 
cons  icier td . 

The  inversions  of  waveforms  for  four  Pahute  Mt-sa  c-xpl os  i  or.s  v,  :■ 
presented  ar.d  discussed  in  some  detail  It:  general,  the  synthetic  v-.vefo;— 
fro::  the  inversions  provided  improvements  over  the  forward  modeling  results  it. 
amplitude,  arrival  time  (except  SCOTCH),  and  in  some  cases,  in  the  detail*  <i 
shapes  of  the  waveforms.  For  BOXCAR,  the  structure  model  obtained  was  net 
significantly  different  than  that  of  Hartzeil,  et  al  .  <  198  3)  .  This  is  not 
particularly  surprising,  since  these  data  were  modeled  in  the  development  of 
the  Hartzeil,  et  al.  (1983)  structure.  For  the  other  events  modeled,  there 
was  verv  little  resolution  of  structure  parameters  deeper  than  about  3  km. 
Above  this,  however,  the  velocity  gradient  near  INLET  is  very  similar  to  that 
near  BOXCAR .  but  the  gradient  is  much  greater  near  MAST.  The  results  from 
SCOTCH  suggest  that  the  water  table  is  somewhat  deeper  than  the  constrain'd 
value  used,  but  that  the  gradient  below  this  is  similar  to  that  of  BOXCAR 
Referring  to  Figure  1,  we  see  that  BOXCAR,  INLET  and  SCOTCH  are  all  located 
well  within  the  Silent  Canyon  Caldera,  the  structural  feature  that  defines 
Pahute  Mesa.  However.  MAST  is  located  on  the  northeast  edge  of  the  mess 
where  one  might  expect  variations  in  velocity  structure. 


The  parameters  of  the  effective  source  functions  were  used  to  develop 
sealing  relations  and  to  estimate  yields  for  the  events  studied.  Barker,  ft 
al .  i  1983a  estimated  vie  Ids  from  the  forward  modeling  results  These  val u-  • 
are  given  in  Table  U.  We  have  obtained  two  estimates  of  vield  from  :!■ 
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:■  vi-rsion  results,  using  relations  determined  for  the  announced  yields  of 
<• : *  her  BOXCAR  or  SCOTCH.  We  prefer  the  estimates  from  BOXCAR ,  which  give 
yields  of  -3’  k:  for  MAST .  391  kt  for  INLET  and  175  kt  for  SCOTCH. 

While  the  resolution  of  the  inversion  results  depends  on  the  range 
distribution  of  the  observations,  and  the  results  are  limited  by  noise  or  the 
•it-  ts  of  lateral  structural  variations  in  the  observed  data,  the  application 
l.if.ute  Mesa  waveforms  has  demonstrated  the  usefulness  of  the  method 
7r  iie-offs  bet  ween  parameters,  particularly  between  structure  and  source 
parameters,  are  quantified  by  the  error  analysis  inherent  in  the  inversion 
method.  Finally,  the  inversion  mav  be  applied  to  near-field  waveforms  from 
ar.y  site,  without  any  a  priori  information  on  the  crustal  velocity  structure. 
This  makes  such  an  inversion  valuable  as  one  of  many  on-site  verification, 
techniques,  enabling  the  calibration  of  test  sites  for  treaty  monitoring 
purposes . 
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EVIDENCE  FOR  LINEAR  MATERIAL  RESPONSE 


AT  HIGH  STRAIN  IN  THE  NEAR  FIELD 


tv* 3W 
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INTRODUCTION 


Some  of  the  most  successful  applications  of  synthetic  seismogram 
techniques  have  occurred  in  modeling  studies  of  near  field  records  from 
explosions  and  earthquakes.  Evidently,  the  crust  in  manv  situations  can  be 
reasonably  approximated  by  a  layered  elastic  model  containing  no  more  than  ten. 
to  twenty  layers.  As  long  as  sources  of  negligible  spatial  extent  are 
considered,  observed  records  are  generally  simple  and  can  be  matched  in  close 
detail  from  the  body  waves  through  the  surface  waves.  Some  examples  of  such 
studies  are  those  of  Helmberger  and  Hadley  (1981),  Burdick  et  al .  (19S4)  and 

Stump  and  Johnson  (1984),  who  modeled  records  from  nuclear  explosions  and 
those  of  Liu  and  Helmberger  (1985),  Hartzell  and  Brune  (1979)  and  Langston  and 
Franco-Spera  (1985),  who  studied  small  earthquakes.  Because  many  crucial 
properties  of  seismic  sources  can  be  constrained  from  such  studies,  there  is  a 
continuing  effort  to  improve  their  resolving  power.  In  order  to  increase  the 
accuracy  of  near-field  source-crust  models  further,  it  is  becoming  necessary 
to  consider  propagation  in  linear  anelastic  media  rather  than  just  linear 
elastic  media. 

Unfortunately,  developing  accurate  models  for  the  distribution  and 
frequency  dependence  of  Q  in  the  shallow  crust  is  just  as  difficult  as 
determining  accurate  Q  models  for  the  aesthenosphere  and  mantle.  Manv 
different  types  of  physical  phenomena  are  believed  to  contribute  a  component 
to  near  field  absorption  including  scattering  (Menke ,  1984),  bulk  losses  due 
to  absorbed  water  (Tittman  et  al.,  1972)  and  losses  due  to  frictional  sliding 
or.  cracks  at  high  strain  (Stewart  et  al . ,  1983;  Mavko ,  1979;  Dav  and  Minster 
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l:-?v'>.  The  last  of  these  presents  a  problem  unique  to  interpretation  of  near 
field  data  because  it  is  only  in  such  data  that  shear  strains  are  large  enough 
to  induce  the  effect.  Figure  17,  which  is  redrawn  from  Mavko  (1979), 
sur.tmarir.es  the  laboratory  data  which  indicates  that  absorption  begins  to 
increase  significantly  in  some  geologic  materials  when  shear  strain  exceeds 
1 0 ’ s .  Two  points  that  are  illustrated  by  Figure  17  which  will  be  significar* 
in  the  the  following  discussion  are  that,  in  general,  the  increase  in 
absorption  is  more  pronounced  in  softer  material,  and  that  the  effect  should 
become  very  apparent  in  such  materials  as  shear  strain  ranges  from  1  to  100 
microstrain. 

In  a  previous  report  (Burdick  and  Barker,  1986),  a  method  was  presented 
for  estimating  shear  strain  fields  from  observed  near-field  velocity  records. 
It  was  utilized  to  measure  shear  strain  levels  for  a  set  of  observations  from 
nuclear  explosions  at  Pahute  Mesa  and  from  a  small  earthquake  in  Imperial 
Valley.  California.  Here  we  present  an  extended  analysis  of  that  data  set  in 
which  we  have  made  a  more  thorough  attempt  to  detect  a  reduction  in  apparent  0 
due  to  high  shear  strain  levels.  We  begin  with  a  brief  review  of  the  methods 
developed  in  the  preceding  study.  Then  we  discuss  the  various  types  of 
evidence  which  might  indicate  reduction  of  apparent  Q  due  to  high  strain  in 
the  Pahute  Mesa  data  set.  This  data  was  recorded  in  a  region  with  a 
dominantly  sandstone  lithology.  Figure  17  shows  that  reduction  of  apparent  Q 
by  high  strain  should  be  most  easily  observed  in  just  such  an  area.  An 
improved  method  for  measuring  apparent  Q  is  presented  and  utilized  to  measure 
the  observed  differences  in  0  between  various  combinations  of  P  and  S  waves 
Correlations  are  sought  between  the  observed  apparent  Q’s  and  the  peak  shear 
strains  induced  by  the  P  and  S  pulses. 
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Figure  17.  Linear  plot  of  attenuation  versus  strain  amplitude  for  six 
different  dry  rocks  (redrawn  from  Mavko  (1979)). 


PREVIOUS  WORK 


In  a  whole  space,  the  strain  pulse  due  to  a  point  source  is  just  the 
vilocitv  pulse  multiplied  bv  the  medium's  slowness.  In  a  lavered  half  space. 

:  hr  velocitv  and  strain  time  histories  differ  to  some  degree.  The  extent  of 
the  difference  can  be  determined  through  synthetic  seismogram  computation. 

The  i.orzoro  components  of  the  stress  tensor  are  computed  automatically  in 
standard  wavenumber  integration  calculations.  Unfortunate lv ,  the  technique 
used  most  often  in  modeling  studies  of  near  field  data  is  summation  of 
generalized  rays  because  it  provides  high  frequency  synthetics  in  much  less 
time.  Some  slight  modifications  of  the  approach  are  required  to  generate 
synthetic  strain-time  histories  rather  than  motion-time  histories.  Following 
Helmberger  (1974),  we  can  express  a  single  generalized  ray  at  high  frequency 
in  the  far  field  as 


ray!,)-'2  ■■  [  l  *  3mf  U  P  ^  ]‘l  1  h) 

\  nr  dt  I  \  q  dt  1  1  ; 


evaluated  along  the  Cagniard  contour  in  the  complex  p  plane.  The  factor  “  i: 
the  product  of  the  reflection  and  transmission  coefficients  along  the  raypath 
A  is  the  source  radiation  pattern  term  and  R  is  called  the  receiver  function 
The  onlv  difference  between  a  strain  generalized  ray  and  a  seismic  motion 
generalized  ray  is  a  change  in  the  expression  for  the  receiver  function. 

These  functions  for  either  motion  or  strain  are  given  in  Table  11.  The 
factors  p  and  7  are  horizontal  and  vertical  slownesses,  a  and  b  are  the  P  and 
S  velocities  and  •  is  plus  or  minus  one  depending  on  whether  the  ray  is  up  or 
downgoing  at  the  receiver  The  functions  on  the  left  are  for  receiver  points 
in  the  medium  and  those  on  the  right  are  for  points  on  the  free  surface. 


Receiver  Functions 


Seismic  Motions 

At  the  Free  Surface 

RFz-injy,- p7)/P7R(p) 

with  R(p)  «  [n,>  ~  PJ)  4^p:f). 

^52  ■  4pr70r7#/^2jR(p) 

^^  =  2r?#(^2-pJ)/^2/?(p) 


The  above  formalism  was  used  to  compute  synthetic  velocity  ar.d  strain 
pulses  for  several  relevant  test  cases.  Definite  differences  existed  in  the 
velocitv  and  strain  time  histories,  but  they  were  small  in  most  cases.  This 
suggested  that  a  method  analogous  to  the  one  that  works  in  a  whole  space  could 
be  devised  for  the  layered  half  space  case;  a  procedure  that  transforms 
observed  velocity  records  to  estimated  strain  records.  Whereas  in  the  former 
case  velocity  is  multiplied  by  slowness  to  obtain  strain,  in  the  latter 
velocity  can  be  transformed  to  strain  by  convolving  the  velocity  pulse  with 
frequency  dependent  transfer  operators.  These  operators  can  be  generated  by 
deconvolving  synthetic  velocity  from  synthetic  strain  records  of  any  given 
source - rece iver  geometry.  In  most  cases,  these  operators  turn  out  to  be  very 
delta-like  functions  with  amplitude  controlled  by  the  velocity  at  the  receiver 
site,  just  as  one  might  expect.  The  first  data  base  processed  using  the 
transfer  operators  was  a  suite  of  near-field  velocity  records  from  6  Pahute 
Mesa  nuclear  explosions.  They  were  SCOTCH,  INLET,  MAST,  HALFBEAK,  ALMENDRO 
and  BOXCAR.  The  recording  instruments  were  L7  velocity  meters,  so  the  signals 
recorded  on  them  were  essentially  velocity  versus  time.  In  order  to  compute 
transfer  functions,  it  was  necessary  to  select  a  layered  velocity  model  for 
Pahute  Mesa.  The  one  chosen  was  presented  by  Barker  et  al .  (1985)  who 
developed  it  specifically  to  give  accurate  near  field  synthetics.  An  example 
of  the  computation  of  a  transfer  operator  for  a  1000  kt  event  located  at  1  km 
depth  in  the  Pahute  Mesa  structure  is  shown  in  Figure  18.  The  synthetic 
vertical  and  radial  velocity  traces  are  shown  on  the  left,  the  four  nonzero 
partial  derivatives  of  displacement  with  respect  to  spatial  coordinates  in  the 
center  and  the  transfer  functions  on  the  right  The  top  strain  is  Ezz,  the 
bottom  is  Epj;  and  the  average  of  the  center  two  is  E^ .  Note  that  the  center 
two  traces  sum  to  zero,  so  the  free  surface  condition  is  maintained.  The 
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Figure  18.  An  example  of  the  computation  of  the  velocity  to  strain  transfer 
functions  for  the  Pahute  Mesa  structure ,  The  velocity  traces  on  the  left  are 
deconvolved  from  the  strain  traces  in  the  center  to  produce  the  transfer 
operators  on  the  right. 
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transfer  functions  are  coirputed  bv  transform!:: m.  the  velocity  pulse  and  the 
corresponding  partial  derivative  pulse  into  the  frequency  domain  using  a  f  a * 
four ier  transform  algorithm  The  latter  is  divided  by  the  former  and  the 
inverse  transform,  taken.  The  signal  to  noise  ratio  of  the  impulses  in  Figure 
IS  is  about  5  to  1 .  To  suppress  the  noise,  a  gauss i an  filter  with  a  cutoff  of 
^  hr  is  applied  to  the  transfer  operator  before  it  is  applied  to  the  data.  Ir 
principle,  both  the  vertical  and  radial  records  could  be  used  along  with  the 
appropriate  transfer  operators  to  estimate  strain.  In  practise,  it  is 
difficult  to  do  so  because  the  synthetics  are  never  really  exact,  and 
generally  the  vertical  is  predicted  better  than  the  radial.  It  is  difficult 
to  achieve  the  delicate  cancellations  required  to  maintain  the  free  surface 
condition.  Therefore,  the  vertical  component  is  used  to  generate  one  of  the 
two  nonzero  strain  components  for  an  explosive  source,  and  the  free  surface 
condition  is  used  to  generate  the  other  from  it. 

The  next  records  to  which  the  velocity  to  strain  transfer  operators  were 
applied  were  from  an  aftershock  of  the  15  October  1979  Imperial  Valley 
earthquake.  The  aftershock  was  studied  in  some  detail  by  Liu  and  Helmberger 
i 1985)  who  provided  a  mechanism,  moment  and  time  function  for  it.  They 
reported  the  event  depth  as  9.5  km.  The  predominant  arrival  at  all  stations 
was  a  strong,  unusually  simple  S  wave,  so  in  this  case  we  computed  S  rather 
than  P  operators.  The  polarity  information  in  the  data  set  was  very  clear, 
and  it  allowed  Liu  and  Helmberger  (1985)  to  constrain  the  mechanism,  to  be 
alm.o st  pure  vertical  strike  slip.  They  reported  a  moment  of  1.0  x  102“ 
dvno-c.Ti  and  a  triangular  time  function  with  a  rise  of  0.1  and  a  fall  of  0.1  s 
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From.  our  modeling  studies,  however,  we  conclude  that  a  moment  value  of  0.6  >. 


10-'  dyne-cm  and  a  time  function  with,  a  rise  of  0.3  and  a  fall  of  0.1  s  are- 
more  accurate . 

The  nonzero  strains  generated  bv  incident  SV  are  E-,  and  .  As  before, 
we  wished  to  insure  that  the  free  surface  condition  would  be  satisfied,  so  we 
used  only  the  radial  component  of  motion  and  generated  both  strains  from  it. 
The  transverse  records  were  used  to  generate  the  third  nonzero  component  of 
strain.  Another  important  difference  between  the  Pahute  Mesa  and  Imperial 
Valley  data  sets  is  that,  in  the  latter,  Liu  and  Helmberger  (1983)  found 
evidence  for  a  very  low  shear  Q  in  the  form  of  a  strong  frequency  shift 
between  P  and  S  waves.  The  values  they  used  for  shear  Q  in  the  top  half  km  of 
their  crustal  model  were  of  the  order  of  10.  This  resulted  in  a  relatively 
distance  independent  t*  of  .132  s.  Values  of  shear  Q  this  low  are  atypical 
and  could  easily  be  interpreted  as  an  indication  that  nonlinear  reduction  of  Q 
due  to  high  strain  is  taking  place.  The  peak  shear  strain  values  for  16 
observing  stations  ranged  from  5  to  65  x  10"6.  We  here  report  on  our 
additional  attempts  to  detect  in  the  Imperial  Valley  and  Pahute  Mesa  data 
sets,  any  evidence  of  nonlinear  material  response. 

EXPLOSIONS 


Devising  means  for  detecting  nonlinear  reduction  in  apparent  Q  from 
near-field  or  teleseism.ic  observations  of  nuclear  explosions  is  difficult 
because  there  are  always  tradeoffs  between  the  effects  of  attenuation  and 
other  processes.  In  the  near  field,  the  decay  rate  of  wave  fields  is 
dependent  on  both  the  elastic  and  the  anelastic  crustal  structure.  At 
teleseism.ic  distances,  estimated  values  for  total  attenuation  (embodied  in  a 
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'ountea  tor.  -n  tr.is  sectior..  we  cc:r;  re  crus:  .1  * 


5r.ru:':ures  and  t*  value:,  estimated  from  nigh  strain  explosion  data  to 
structures  and  t*s  estimated  bv  other  means. 


In  the  near  field  modeling  studies  of  Burdick  et  al  .  (,  198-1;  and  Hartz 


al  . 1  1983}  .  elastic  crustal  stn 


developed  for  Amchitka  and 


i  Cd.ute  Mesa  by  forward  modeling  of  observed  velocity  pulse  amplitudes  and  v 
shapes.  The  question  is  whether  nonlinear  effects  v:ere  inadvertently  mapped 
into  details  of  the  crustal  model.  Figure  19  shows  the  estimated  peak  shea 
strain  for  the  Pahute  Mesa  data  as  a  function  of  range.  Theoretical 
strain-  distance  curves  computed  assuming  the  iiartzell  et  al  .  (  1933.)  structu 
and  elastic  material  response  for  the  largest  (BOXCAR,  1300  kt)  and  smalles 


SCOTCH.  135  kt)  events  are  also  shown .  The  theoretical  strain  for  the 
smallest  event  does  not  fall  below  one  micros' rain  until  25  km .  Thus,  the 
on.'  ire  Pahute  Mesa  near  -  field  data  set  is  within  the  strain  regime  in  which 
tht  laboratory  data  indicates  that  nonlinear  processes  might  be  significant 

The  Amchitka  event,  MILROVJ,  was  comparable  to  BOXCAR  in  size  so  it  probabl v 

.v-re rat ed  shear  strains  of  the  same  order.  The  other  Amchitka  event, 

I  y V.  was  three  or  I  <  r  times  as  large,  but  the  velocity  pulses  it 

need  were  very  r.earl  v  the  same-  size  as  tl  ->r-e  from  MILROW .  The  two  eve 

'  ■:>)  :  •  probably  gerv  ia'cd  comparable  shea:  '.trains  at  the  free  surface 
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igure  19.  The  decay  of  peak  shear  strain  near  to  nuclear  explosions.  The 
argest  observed  value  of  shear  s  rain  is  10"3.  Theoretical  curves  are  shown 
or  the  largest  and  smallest  events  studied.  The  observations  appear 
olio w  the  theoretical  predictions 
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the  vears  .  both  I'S  n.uc  I  ear  test  sites  have  beer,  character ;  zed  ir, 
.'hvsieal  de:  ail  Ve  1  or  i  t  v  profile;,  b.i.sed  on  low  st  rain  tv  pcs  of 
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crust  moce.s 


of  the  previous  section  There  is  certainly  no  evider. 


to  ftgr.fs:  that  ar.v  biasing  has  occurred.  The  observed  and  predicted 
amplitude  deeav  with  distance  is  given  in  figures  9,  11  arid  13  for  both  the 
forward  and  inverted  crus:  models  There  is  i.o  indication  that  the  wave  f  i t- 1  i 
is  being,  anomalously  attenuated  with  distance.  Crustal  structure  models  for 
Am.chi  :k.«  based  or.  hir’d,  s  t  r i  r.  versus  low  strain  data  are  compared  in  Figure 
1 0  The  MILR'w  and  LANN  I F. ;  N  node’.':  are  based  on  near  field  observations,  arid 
t  he  Eng  dab.  1  1  ■5  "c>  model  is  based  on  low  strain  refraction  information. 

Again,  there  appears  tc  be  no  bias.  Figures  21  and  22  show  the  observed 
versus  predicted  peak  velocity  value  as  a  function  of  range  for  MILROV  and 
CANNIKIN  respectively.  There  is  no  evidence  for  a  high  rate  of  attenuation 
that  is  not  accounted  for  in  the  linear  elastic  model.  Note  that,  for  the 
vertical  component,  the  theoretical  curve  continues  to  predict  the 
observations  even  inside  the  spall  zone  (less  than  3  km)  though  this  is  not 
true  for  the  radial.  In  any  event,  there  is  no  evidence  in  the  comparisons  of 
high  strain  versus  low  strain  crustal  models  that  any  biasing  due  tc  high 
strain,  effects  has  occurred.  Likewise,  there  is  no  evidence  of  low  apparent  C 
in  the  peak  ampl i tude-distance  data 


It  is  difficult  to  assess  how  a  substantial  unaccounted  for  reduction  in 
effective  0  in.  the  near  field  would  affect  at.  estimate  of  t*  based  on 


t.  e  1  “  st-  i  sir  i  c  amp  1  i  t  udes 


isual  assumption  in.  such  measurements  is  that  the 


explosive  source  is  effect iv<  ly  spherical  If  this  assumption  remains  true, 
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Figure  21.  A  comparison  of  the  observed  velocity  data  and  elastic  rr.odel 
predictions  for  MILROV.  Observed  peak  vertical  (circles)  and  radial  (squares) 
velocities  and  elastic  model  predictions  (solid  and  dashed  lines)  plotted  as  a 
function  of  epicentral  distance. 


bur  the  reduction  in  Q  extends  bevond  the  near  field  data,  the  initial 
strength  of  the  far  field  elastic  wave  field  as  measured  from  that  data  would 
be  overestimated.  The  average  value  of  t*  for  the  whole  ravpath  would  be 
overestimated.  If  the  tone  in  which  effective  Q  is  reduced  is  very  aspheric  : 
and  confined  to  the  surface  1  avers ,  the  initial  field  strength  would  be 
underestimated  and  t*  for  the  whole  path  underestimated  The  only  effect 
which  should  be  present  in  either  case  would  be  a  yield  (peak  shear  strait.) 
dependence  of  the  phenomenon. 

As  noted  above,  the  important  tradeoff  in  such  measurements  is  between 
uncertainties  in  the  yield  scaling  of  the  RDP  and  t* .  The  apparent  yield 
scaling  for  the  Arrthitka  test  site  based  on  near*field  modeling  was  given  bv 
Burdick  et  al .  (1984),  and  for  Pahute  Mesa  bv  Burger  et  al .  (1987).  In  both 

cases,  the  RDP  formalism  used  was  the  one  proposed  by  Helmberger  and  Hadlev 
,1981).  Another  widely  recognized  RDP  and  yield  scaling  model  for  Pahute  Mes 
is  the  one  presented  by  Mueller  and  Murphy  (197]).  This  model  was  also  based 
to  some  extent  on  near  field  field  data  but  primarily  on  near  regional  data 
where  strains  are  much  lower.  The  RDP  spectra  are  compared  in  Figure  23  for 
large  (BOXCAR!,  intermediate  (HALFBEAK.)  and  small  (SCOTCH)  event.  The 
differences  at  long  period  are  not  significant  since  the  near  field  data 
contained  no  long  period  information.  The  corner  frequencies  and  spectral 
falloff  rates  near  the  peak  track  very  well  over  the  entire  range  of  yields 
As  explained  by  Burger  et  al .  (1987),  the  absolute  offset  of  the  two  spectra 

results  in  the  Mue 1 le r - Murphy  source  giving  consistently  higher  t*  estimates, 
hut  that  is  not  important  to  this  discussion  The  point  is  that  the  yield 
staling  of  the  RDP  indicated  by  the  near  field  data  does  not  appear  to  be 
biased  in  anv  significant  way.  The  average  earth  t*  values  determined  by 
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Figure  23.  Predicted  far-field  source  displacement  spectra  for  BOXCAR . 
HALFBLAK,  and  SCOTCH.  The  Mueller-Murphy  (1971)  source  models  are  based  on 
the  announced  yields.  The  Helmberger-Hadley  (1981)  source  models  are  based  on 
the  near-field  modeling  results  of  Barker  et  al,  (1985). 
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or  of  four  or  five. 


Pab.ute  Mesa,  rhere  is  a  total  variation  of  about  30*  in  estimated  t ‘  with 


variation  in  yield  of  an  order  of  magnitude  2 he  extremal  values  were  for 


largest  and  smallest  events  ,  but  other  than  that  there  are  no  convincing 


svstetatics .  There  is  certainly  no  strong  evidence  for  reduction  of  effective 


Q  in  the  near  field  having  a  large  effect. 


EARTHQUAKES 


"he  aftershock  of  the  Imperial  Earthquake  occurred  much  deeper  than  the 


typical  depth  of  explosions.  In  fact,  it  was  so  deep  that  the  body  wave 


propagation  to  the  near  field  stations  remained  very  simple.  The  arrivals  are 


well  represented  by  a  single  upgoing  generalized  ray.  The  rays  to  all 


stations  depart  nearly  horizontally  from  the  source  and  then  turn  sharply 


upward  under  the  stations  because  of  the  strong  velocity  gradient.  The 


observed  S  wave  pulses  are  very  simple  in  character,  and  easy  to  model  as 


shown  in  Figure  7U  .  The  SH  and  SV  waveshapes  are  very  similar  because  the 


rays  emerge  nearly  vertically  at  the  receiver.  The  absence  of  ringing 


indicates  that  reverberations  at  the  stations  are  not  significant.  As  noted 


above,  the  peak  shear  strains  generated  by  the  16  S  waves  shown  were  estimated 


in  previous  work,  and  our  purpose  here  is  to  determine  whether  there  is 


significant  evidence  of  reduction  of  effective  Q  due  to  these  strains. 


The  observed  decay  of  peak  shear  with  distance  for  the  aftershock  is 


shown,  in  Figure  25.  The  theoretical  curves  for  a  large  and  small  explosion 


! torn  Figure  19  are  shown  for  reference.  The  levels  of  shear  strain  generated 
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Table  12  -  Measured  Values  of  t*  versus  Yield 


Test  Site 

Event 

Yield 

t* 

(kt) 

(sec 

AMCHITKA 

CANNIKIN 

<  5000 

.94 

MILROW 

1000 

.85 

NTS 

BOXCAR 

1300 

.92 

ALMENDRO 

570 

.69 

MAST 

520 

.80 

INLET 

500 

.72 

HALFBEAK 

365 

.  74 

SCOTCH 

155 

.61 
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bv  the  aftershock  are  quite  comparable  to  those  generated  bv  the  nuclear 
events  although  the  shear  modulus  of  the  Imperial  Valiev  sediments  is  lower 
than  those  at  the  surface  of  Pahute  Mesa.  The  decay  of  strain  with  distance 
is  different  for  the  explosion  and  earthquake  cases,  but  this  is  primarily  du-. 
to  the  greater  depth  of  the  aftershock.  There  is  more  scatter  in  the 
earthquake  data,  but  this  is  at  least  in  part  caused  by  the  azimuthal 
radiation  pattern  of  the  double  couple  which  has  not  been  corrected  for  in 
this  Figure.  Given  the  evidence  for  the  high  peak  shear  strains  in  Imperial 
Valley  and  the  susceptibility  of  sandstone  lithologies  to  reduction  in 
apparent  Q  by  strain,  the  Imperial  Valley  aftershock  data  set  appears  to  be 
ideal  for  searching  for  evidence  of  this  phenomenon. 

Liu  and  Helmberger  (1985)  initially  concluded  that  the  effective  Q  in  the 
Imperial  Valley  had  to  be  relatively  low  because  of  the  strong  shift  in 
frequency  content  between  P  and  S  waves.  Such  a  shift  has,  of  course,  several 
possible  origins.  These  include  elastic  propagation  differences  for  P  versus 
S  waves  and  directivity  differences  at  the  source  as  well  as  strong 
differential  attenuation.  The  elastic  propagation  is  so  simple  in  its 
geometry,  however,  that  it  is  not  a  likely  candidate.  Source  directivity 
effects  can  be  damped  by  averaging  over  azimuth,  so  it  should  be  possible  to 
isolate  the  effects  of  differential  attenuation  to  some  degree.  If  it 
correlates  with  the  difference  in  peak  shear  strain  carried  by  P  and  S,  we 
would  have  some  evidence  of  nonlinear  material  response.  Figure  26  shows  some 
examples  of  the  strong  P  to  S  frequency  shift.  Three  examples  of  observed  P 
waves  are  shown  in  the  left  column.  Even  though  they  are  very  high  in 
frequency  content,  they  are  relatively  simple  for  near-field  S  waves  just  as 
the  S  waves  in  Figure  25  are.  In  the  center  column  are  Futterman  attenuation 
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Figure  26.  Trans f orma t i on  of  P  into  S- waves  wit  a  4 
The  observed  P- waves  are  shown  on  the  lei:  the  C  ‘ 

center  column  and  the  filtered  P  is  compared  to  • 

the  variability  of  the  apparent  relative  a : • en  ;  . c- 
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operators  that  transfer  the  P  pulses  into  pulses  that  most  resemble  the  S 
waves.  The  criterion  for  selecting  these  operators  is  discussed  in  the 
following  section.  The  filtered  P  is  compared  to  the  associated  S  pulse  in 
the  right  column.  The  peak  shear  strain  carried  by  each  pulse  is  also 
indicated.  The  amount  of  relative  attenuation  is  quite  variable,  and  in  these 
examples  it  appears  as  if  it  may  be  related  to  the  difference  in  P  and  S 
strain.  Also,  the  filtered  P  and  S  waves  match  each  other  in  close  detail 
indicating  that  a  difference  in  attenuation  might  be  the  most  significant 
difference  between  them. 


It  must  be  acknowledged  at  this  point  that  representing  the  effects  of 
nonlinear  material  response  with  simple  attenuation  operator  may  turn  out  to 


o 


be  a  very  poor  approximation.  As  discussed  in  day  and  Minster  (1986),  it  is 
necessary  to  consider  propagation  in  media  with  generalized  constituative  laws 
to  treat  the  problem  in  its  entirety.  However,  solutions  to  such  problems  are 
much  more  complex  than  those  commonly  used  in  seismology,  and  it  is 
appropriate  to  first  search  for  the  types  of  correlations  in  observations  that 
might  be  expected  if  nonlinear  effects  are  important.  If  such  patterns  begin 
to  emerge,  they  can  be  modeled  more  exactly  at  a  later  time. 


An  I-proved  Me thod  of  Measuring  t* 


Measuring  the  value  of  t*  from  the  relative  frequency  content  of  P  and  S 
waves  is  a  commonly  used  technique.  Kanamori  (1967),  who  used  the  method  to 


measure  mantle  t*  from  ScP  and  ScS  phases,  discusses  the  theory  and 
assumptions  behind  the  approach.  For  reasons  of  mathematical  convenience,  ht 
suggests  choosing  an  optimal  value  for  t*  by  taking  the  natural  log  of  the 


spectral  ratio  of  S  to  P  and  fitting  a  straight  line  to  the  results.  In  this 
discussion,  we  will  refer  to  this  as  measuring  t*  by  minimizing  the  Ln  Least 


Squares  or  LLS  norm.  That  is, 


„  r 


LLS  =  min 


In (5/?) - 


co  t 


d  co 


(30) 


where  S  and  P  are  the  amplitude  spectra.  Using  the  definition  for  a  frequency 
independent  attenuation  operator  A, 


A(t*)  =  exp 


(30) 


the  LLS  norm  can  be  recognized  as 


LLS  =  m  in 


[ln(S)-lnM(f‘)P)]2dco 


(31  ) 


Note  that  because  only  the  amplitude  spectra  are  used  the  phase  information  is 
neglected . 

An  alternate  norm  that  is  frequently  used  to  analytically  compare  pulses 
is  the  one  suggested  by  Burdick  and  Mellman  (1974)  .  Because  they  express 
their  norm  in  terms  of  cross  correlation  operators,  many  fail  to  realize  that 
'  it  is  just  a  standard  least  squares  (LS)  norm.  It  could  be  expressed  in  this 

I 

' 

case  as 

/S-minJlS-.1(f*)/’l?rft  (30  j 
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where  the  S  and  P  are  aligned  in  time  at  max  correlation  and  where  each  pulse 
is  normalized  to  unit  squared  area.  The  squared  difference  function  in  the 
brackets  is  just  a  function  of  time,  so  an  immediate  application  of  Parseval's 
theorem  to  it  gives 

LS  =  min  j  [S  -  A(t*)P}7duj  (33) 

where  obviously  P,  S  and  A  are  to  be  in  the  time  domain  in  the  first  of  these 
two  equations  and  in  the  frequency  domain  in  the  second.  Thus  the  LLS  and  the 
LS  norms  differ  only  in  that  one  minimizes  the  difference  between  natural  logs 
of  spectra  and  the  other  the  spectra  themselves,  and  that  the  former  ignores 
phase  information  while  the  other  utilizes  it.  While  the  LLS  norm  is  easier 
to  apply,  it  does  have  some  undesirable  numerical  properties  as  we  show  in  the 
following . 

To  illustrate  the  difference  in  the  two  norms,  we  consider  a  sample 
record  from  the  Imperial  Valley  data  set  (ELCE  Radial  S  and  Vertical  P) .  The 
LLS  norm  indicates  that  t*  is  .03  s  and  the  LS  norm  indicates  .09  s.  The  LLS 
line  was  fit  over  the  frequency  range  of  1  to  10  hz  since  this  is  the  range 
where  signal  to  background  noise  is  larger  than  1.  The  fit  of  the  LLS  line  to 
the  ln(spectral  ratio)  is  compared  to  the  fit  of  the  LS  line  in  Figure  27. 

The  LS  line  has  been  arbitrarily  given  the  same  value  at  1  hz  as  the  LLS  line. 
It  is  clear  that  the  LLS  norm  fits  the  5  to  10  hz  information  better  than  the 
LS  norm  presumably  because  it  is  more  sensitive  to  it.  However,  one  would 
wish  for  a  norm  to  be  most  sensitive  to  the  data  that  is  most  reliable,  and  ir. 
general  that  is  not  the  highest  frequency  information.  To  further  establish 
the  strong  sensitivity  of  the  LLS  norm  to  the  higher  frequencies,  it  was 


applied  to  a  suite  of  12  Imperial  Valley  aftershock  record  pairs;  first 
fitting  1  to  10  hz  and  then  1  to  15  hz  data.  On  average,  including  the  extra 
high  frequency  information  changed  the  estimated  value  of  t*  by  a  factor  of  3 
The  estimated  value  for  the  spectral  slope  is  very  unstable  and  controlled  b\ 
high  frequencies.  The  average  estimated  value  was  about  the  same  in  the  two 
cases  since  as  many  slopes  were  perturbed  up  as  down.  Spectra  and  spectral 
ratios  are  always  noisy,  so  many  have  become  accustomed  to  poor  fits  between 
theory  and  data,  but  it  is  worthwhile  to  note  that  the  spectral  ratio  in 
Figure  27  does  not  really  look  very  much  like  a  straight  line 


Figure  28  compares  the  fits  of  the  two  norms  in  the  time  and  linear 
spectral  domains.  The  S,  filtered  P  and  difference  functions  are  shown  in 
time  and  the  spectra  of  the  two  difference  functions  are  plotted  with 
frequency.  The  solid  curve  (t*-.09  s)  clearly  has  less  area  under  it.  The  LS 
norm  minimizes  the  squared  area  under  the  difference  function  curves  in  either 
domain.  The  norm  is  most  sensitive  to  those  frequencies  with  largest  spectral 
amplitude.  Unless  the  fit  between  the  two  time  functions  is  extraordinary, 
the  difference  function  generally  has  high  amplitude  at  the  same  frequencies 
the  signals  do.  The  largest  amplitudes  in  this  case  are  between  1  and  5  hz . 

It  can  be  seen  in  the  time  domain  that  this  is  where  most  of  the  signal  energv 
is.  At  frequencies  higher  than  this,  the  norm  only  requires  that  the  spectra 
be  small.  Considering  the  logarithms  of  functions  takes  information  down  in 
amplitude  by  orders  of  magnitude  and  gives  it  equal  or  in  this  instance  even 
more  weight  than  higher  amplitude  information.  Based  on  these  tests  we 
conclude  that  the  LS  norm  is  superior  for  estimating  t* . 
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Figure  28  A  comparison  of  the  difference  functions  obtained  by  the  least 
squares  norm  (,09s)  versus  that  prescribed  by  the  Ln  least  squares  norm 
(,09s).  The  top  time  trace  is  observed  radial  S  at  ELCE,  the  next  is  filtered 
P  and  the  third  is  the  time  domain  difference  function.  The  solid  line  is  the 
spectrum  of  the  difference  function  for  least  squares  and  the  dashed  line  for 
Ln  least  squares. 


The  Relation  Be  tween  Relative  Strain  and  Attenuation  at  Imperial  Valiev 

As  discussed  above,  there  are  two  strong  difficulties  with  trving  to 
detect  nonlinear  reduction  in  Q  by  correlating  relative  attenuation  of  P  and  S 
waves  with  the  relative  shear  strains  that  they  induce.  The  first  is  that 
source  directivity  causes  frequency  shifts  between  P  and  S  which  can  not  be 
easily  distinguished  from  the  effects  of  relative  attenuation.  In  order  to 
suppress  this  effect,  it  is  necessary  to  average  observations  from  many 
azimuths.  The  second  is  that  the  relationship  between  effective  Q  and 
relative  shear  is  likely  to  be  a  very  complex  one.  All  that  can  be  attempted 
is  to  find  correlations  between  the  two  variables  by  considering  large  data 
sets.  The  Imperial  Valley  aftershock  data  set  contains  6  good  P  wave  and  32  S 
wave  observations.  Given  the  need  for  strong  averaging,  we  elected  to  measure 
the  relative  attenuation  of  each  P  wave  with  every  S.  This  yielded  192 
separate  observations.  The  LS  norm  defined  above  yields  and  estimate  of 
relative  t*  rather  than  a  direct  estimate  of  shear  Q.  However,  we  can 
determine  the  latter  because  we  have  a  reliable  crustal  velocity  model  for 
Imperial  Valley.  The  elastic  compressional  and  shear  wave  travel  times  are 
computed  from  this  model  and  shear  Q  is  determined  from 
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The  value  of  c  varies  through  the  crust,  but  we  repeated  our  calculations  for 
the  minimum  and  maximum  values  and  found  that  it  does  not  have  a  significant 
effect 

The  results  of  the  measurements  are  shown,  in  Figure  29.  As  might  be 
expected  in  situations  in  which  there  are  many  uncontrolled  variables,  the 
data  shows  a  large  amount  of  scatter.  On  the  other  hand,  most  tvpes  of  0 
measurements  show  significant  scatter.  It  should  be  noted  that  the  vertical 
scale  is  linear  and  not  logarithmic  and  that  virtually  all  values  are  confined 
between  0.0  and  0.05.  The  trend  of  Q'1  with  shear  strain  measured  in  the 
laboratory  as  shown  in  Figure  17  is  drawn  in  for  reference.  A  regression  line 
which  was  fit  to  the  data  is  also  shown.  It  shows  a  slight  positive  slope, 
but  that  slope  has  a  low  level  of  statistical  significance.  In  other  words, 
there  is  certainly  not  strong  evidence  that  effective  Q  is  being  reduced  due 
to  high  shear  strains  though  there  is  perhaps  some  suggestion  of  it.  It  is 
worth  noting  that  the  value  of  shear  Q'1  at  low  strain  in  Figure  17  is  .012 
which  is  very  near  to  the  average  of  this  suite  of  measurements  which  lends 
some  additional  confidence  as  to  the  validity  of  the  measurement  technique. 

One  final  type  of  information  in  the  earthquake  data  set  which  might  be 
looked  at  to  detect  anomalous  reduction  in  shear  Q  is  the  observed  amplitudes 
of  wave  fields  with  respect  to  the  elastic  predictions.  This  is  somewhat  more 
difficult  for  a  double  couple  than  an  explosion  source  because  of  the 
am ituthal  radiation  pattern  More  distant  points  may  undergo  stronger  shear 
strains  than  nearer  ones  if  the  nearer  ones  are  on  nodes.  The  final  Figure. 
Figure  30.  shows  the  ratio  of  the  observed  peak  shear  strain  over  the  values 
predicted  from  elastic  theory  versus  observed  peak  shear  strain.  Again,  there 
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Figure  29.  The  inverse  of  shear  Q  as  a  function  of  the  relative  shear  strains 
carried  by  the  P  and  S  pulses.  The  higher  slope  straight  line  is  a  fit  to 
laboratory  observations  and  the  lower  is  a  regression  fit  to  the  data. 
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Figure  30.  The  ratio  of  observed  to  predicted  (linear  elastic)  peak  shear 
strain  for  stations  that  recorded  the  Imperial  Valley  aftershock  as  a  function 
of  the  observed  peak  shear  strain.  There  is  no  evidence  that  signals  are 
anomalously  attenuated  as  peak  shear  increases. 


is  a  slight  trend  to  the  data,  but  nothing  strongly  convincing.  Also  the 
sense  of  the  trend  is  opposite  to  what  would  be  expected  if  more  attenuation 


occurs  at  higher  shear  strain. 


DISCUSSION 

There  are  two  questions  which  need  to  be  addressed  with  regard  to 
nonlinear  reduction  in  Q  due  to  high  strains  in  the  seismic  near  field.  The 
first  is  can  it  be  observed  to  occur  at  all,  and  the  second  is,  if  it  does 
occur ,  how  significant  is  it  to  interpretation  of  near  field  data?  Several 
authors  have  argued  on  both  theoretical  and  observational  grounds  that  it 
never  occurs  in  the  fieid  (Savage,  1969;  Winkler  et  al . ,1979)  though  they  may 
not  have  been  considering  motions  near  to  large  nuclear  explosions.  The 
answer  to  the  first  question,  as  far  as  this  study  is  concerned  is  that  any 
evidence  for  it  is  weak.  There  was  the  difference  in  the  SCOTCH  and  BOXCAR  t* 
measurements  in  Table  12  and  the  trend  in  Figure  27.  However,  it  could  be 
argued  that  the  effect  should  be  very  slight  in  data  of  the  various  types 
examined  here.  The  effect  is  strongly  pressure  sensitive  and  may  only  occur 
in  shallow  layers  near  the  surface  of  the  earth.  The  answer  to  the  second 
question  is  that  accounting  for  this  phenomenon  is  almost  certainly 
unimportant  to  the  correct  interpretation  of  near-field  data.  Linear  elastic 
or  linear  anelastic  models  are  all  that  is  required  to  explain  all  important 
observed  phenomena  in  the  range  of  2  to  20  km.  It  is  interesting  to  find  that 
small  earthquakes  generate  strains  comparable  to  those  generated  by  nuclear 
explosions  even  if  those  events  are  as  large  as  1000  kt  in  size.  The 
widespread  success  of  linear  anelastic  models  in  explaining  near  field  data 
from  earthquakes  should  lend  confidence  to  the  interpretation  of  near  field 
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